The use of human liver microsomes as a model system to evaluate the impact of complex genomic traits (i.e., linkage-disequilibrium patterns, coding, and non-coding variation, etc.) on drug metabolism remains challenging. To overcome these challenges, we propose the use of non-linear mixed effect models (NLME) as an unconventional approach to evaluate the impact of complex genomic traits on the metabolism of xenobiotics in vitro. In this in-silico study, we present a practical use case for such an approach using previously published in-vitro CYP2D6 data and explore the impact of sparse sampling, and experimental error on known kinetic parameter of CYP2D6 mediated formation of 4-hydroxy-atomoxetine in human liver microsomes.
This supplement provides the R code used to conduct the
in-silico portion of the study, as well as useful functions for
evaluating/visualizing results.
# R Packages -----
library(nlme) # Non-Linear Mixed Effect Modeling
library(tidyverse) # Data Manipulation and Visualization
library(ggeasy) # Wrapper for Plot Manipulation
library(ggpubr) # Additional Tools for Plot Configuration
library(scales) # Additional Tools for Adjusting Plot Scales
library(cowplot) # Data Visualization Themes and Tools
library(kableExtra) # Additional Tools for Creating Custom Tables
# Package References -------
# knitr::write_bib(c(.packages()), "References/packages.bib")Citation information for R packages1–16 used for data analysis can be found in the References section.
CYP2D6 mean estimates and frequencies (Dinh et. al., 2016) were imported and summarized by genotype (referred to as diplotype in R script). Additionally, an Eta table and parameter estimate table were also generated to be reference later on in the data analysis.
# Atomoxetine Dataset (Dinh et al., 2016)
atx.data <- readxl::read_xlsx("Input Data/Atomoxetine Data - 2016.xlsx")
# Summary of Atomoxetine Data by genotype ----
diplotype.data <- atx.data %>%
group_by(Diplotype) %>%
# Summarize mean estimates of the data
summarise(Vmax = mean(Vmax),
Km = mean(Km),
Score = unique(`Activity Score`),
n = n()) %>%
ungroup() %>%
mutate(Freq = n/sum(n),
Vmax.Eta = ifelse(Diplotype != "1/1",
round(Vmax-Vmax[Diplotype == "1/1"],
digits = 2),
round(Vmax, digits = 2)),
Km.Eta = ifelse(Diplotype != "1/1",
round(Km-Km[Diplotype == "1/1"],
digits = 2),
round(Km, digits = 2))) %>%
select(-n)
# Eta Table (Deviations from the reference group) -----
actual.eta <- diplotype.data %>%
select(Diplotype, Vmax.Eta, Km.Eta) %>%
rename("Vmax" = "Vmax.Eta",
"Km" = "Km.Eta") %>%
pivot_longer(Vmax:Km, names_to = "Parameter",
values_to = "Actual") %>%
arrange(desc(Parameter))
# mean estimates of Vmax and Km for each genotype -----
actual.est <- diplotype.data %>%
select(Diplotype, Vmax, Km) %>%
pivot_longer(Vmax:Km, names_to = "Parameter",
values_to = "Actual") %>%
mutate(Actual = round(Actual, digits = 2)) %>%
arrange(desc(Parameter))A custom simulation function population.sim() was
scripted to generate a virtual population with a pre-defined
distribution of CYP2D6 genotypes, and to simulate
Michaelis-Menten kinetics for each individual. The inputs for the
function the population.sim() include the following:
n - number of subjects per genotype group (default = 10).
CV%_D - Average standard deviation for each genotype group as a percentage of the true estimate for Vmax and Km (between subject variability for each genotype group) (default = 25).
CV%_V - Average coefficient of variation across all points for simulated Michaelis-Menten kinetics for each individual (residual error).
Start - Starting concentration from Michaelis-Menten simulation (default = 0).
End - Ending concentration for Michaelis menten Kinetics (default = 100 uM).
By - Simulated concentration increment (default = 0.5).
pop_freq - TRUE or FALSE argument indicating whether to use true population frequency (default is FALSE).
seed - Random number generator seed (for reproducibility, default = 23457).
data - Dataframe to be used as reference for
population averages. Must include genotype, Vmax, Km, and Frequency
columns (default = diplotype.data).
Within-group between subject variability was set at 25% CV for each genotype (one standard deviation = 0.25 x mean of true value). Below is a simulated population containing 9000 individuals (1000 per diplotype groups) with mean and standard deviations for each group parameter.
population.sim <- function(
# Pre-define number of subjects per diplotype group
n = 10,
`CV%_D` = 25,
`CV%_V` = 0,
# Pre-define Michaelis Menten Kinetic Settings
Start = 0,
End = 100,
By = 0.5,
# Turn off or on use of true population frequency
## If set to true "n" will equal total population number
pop_freq = FALSE,
# Pre-define reference data, and set seed
data = diplotype.data,
seed = 23456){
# ======== Create Data Containers ========
datasets <- list()
pop.data <- tibble()
# ======== Create Datasets by Diplotype ========
for (i in seq_along(data[[1]])){
# Setting Population Specific Frequency
n_pop <- if_else(pop_freq == FALSE, n,
ifelse(pop_freq == TRUE,
round(n*data[[i, c("Freq")]]), 0))
# Specify Diplotype for current loop
Diplotype = rep(data[[i,1]], n_pop)
# Random assignment of Vmax values N(0,sigma)
set.seed(seed)
Vmax_eta <- rnorm(n_pop, sd = `CV%_D`)
Vmax = data[[i,2]] + (Vmax_eta/100)*data[[i,2]]
# Random assignment of Km values N(0,sigma)
set.seed(seed)
Km_eta <- rnorm(n_pop, sd = `CV%_D`)
Km = data[[i,3]] + (Km_eta/100)*data[[i,3]]
# Store each diplotype group in a list container
temp <- tibble(Diplotype, Vmax, Km)
datasets[[i]] <- temp
}
# ======== Combine Diplotype Datasets ========
for (i in seq_along(data[[1]])){
pop.data <- rbind(pop.data, datasets[[i]])
}
# ======== Simulate Michaelis Menten ========
set.seed(seed)
pop.data <- pop.data %>%
mutate(ID = factor(seq_along(Diplotype),
levels = seq_along(Diplotype))) %>%
relocate(ID) %>%
expand_grid(S = seq(Start, End, By)) %>%
mutate(V = round((Vmax*S)/(Km+S), digits = 2))%>%
group_by(ID) %>%
# Additon of Residual Error ------------
mutate(resid_pct = round(rnorm(n = length(S),sd = `CV%_V`),
digits = 3),
V_adj = V+(V*resid_pct/100)) %>%
ungroup()
return(pop.data)
}dist.plot <- ggplot(
data = population.sim(n=1000,`CV%_D` = 25, Start = 0.5,
By = 5, pop_freq = F) %>%
group_by(ID) %>%
filter(S == min(S)),
aes(x = Vmax))+
geom_histogram(aes(fill = Diplotype), color = "black", bins = 100)+
ggeasy::easy_remove_legend_title()+
theme_bw(base_size = 14)+
xlab(bquote(V[Max]))+
ylab("Population (n)")
dist.grb <- dist.plot +
ggeasy::easy_remove_legend() +
scale_x_log10()+
xlab("Log Scale")+
ylab("Count (n)")
dist.plot2 <- ggplot(
data = population.sim(n=1000,`CV%_D` = 25, Start = 0.5,
By = 5, pop_freq = F) %>%
group_by(ID) %>%
filter(S == min(S)),
aes(x = Km))+
geom_histogram(aes(fill = Diplotype), color = "black", bins = 100)+
ggeasy::easy_remove_legend_title()+
theme_bw(base_size = 14)+
xlab(bquote(K[M](uM)))+
ylab("Population (n)")+
scale_x_log10()
dist.grb2 <- dist.plot +
ggeasy::easy_remove_legend() +
scale_x_log10()+
xlab("Log Scale")+
ylab("Count (n)")
## Vmax and Km Distribution Plots
ggpubr::ggarrange(dist.plot +
annotation_custom(ggplotGrob(dist.grb),
xmin = 1000, xmax = 2500,
ymin = 250, ymax = 1250),
dist.plot2, ncol = 1, common.legend = T, legend = "right")Distribution of Michaelis-Menten Parameters for Virtual Population (n=9000) stratified by CYP2D6 genotype. Km represent the Michaelis-Menten constant, and Vmax represents maximal 4-OH Atomoxetine formation (pmol/min/mg protein) by CYP2D6 in human liver microsomes.
population.sim(n = 1000, End = 0,`CV%_D` = 25) %>%
group_by(Diplotype) %>%
summarise(Vmax_mean = round(mean(Vmax),2),
Vmax_sd = round(sd(Vmax),2),
Km_mean = round(mean(Km),2),
Km_sd = round(sd(Km),2),
`n` = length(Vmax),
Vmax = paste0(Vmax_mean," ± ",Vmax_sd),
Km = paste0(Km_mean," ± ",Km_sd)) %>%
ungroup() %>%
left_join(diplotype.data %>%
mutate(across(where(is.numeric), round, 2)) %>%
select(Diplotype, Vmax, Km) %>%
rename("Actual (Vmax)" = "Vmax", "Actual (Km)" = "Km"),
by = "Diplotype") %>%
select(Diplotype,`Actual (Vmax)`, Vmax, `Actual (Km)`,Km, n) %>%
rename("Simulated." = "Vmax",
"Simulated" = "Km",
"Actual." = "Actual (Vmax)",
"Actual" = "Actual (Km)") %>%
kbl(table.attr = "style='width:60%;'",
caption = "Simulated Population Parameter Estimates") %>%
kable_classic("striped", full_width = T) %>%
add_header_above(c(" ", "Vmax Estimate" = 2, "Km Estimate" = 2, " "))|
Vmax Estimate
|
Km Estimate
|
||||
|---|---|---|---|---|---|
| Diplotype | Actual. | Simulated. | Actual | Simulated | n |
| 1/1 | 354.50 | 352.17 ± 88.91 | 1.45 | 1.44 ± 0.36 | 1000 |
| 1/2 | 470.33 | 467.24 ± 117.97 | 1.39 | 1.38 ± 0.35 | 1000 |
| 1/2x2 | 1410.00 | 1400.73 ± 353.65 | 1.82 | 1.81 ± 0.46 | 1000 |
| 1/4 | 274.00 | 272.2 ± 68.72 | 0.92 | 0.91 ± 0.23 | 1000 |
| 2/2 | 252.00 | 250.34 ± 63.2 | 4.59 | 4.56 ± 1.15 | 1000 |
| 2/3 | 251.00 | 249.35 ± 62.95 | 1.32 | 1.31 ± 0.33 | 1000 |
| 2/4 | 95.90 | 95.27 ± 24.05 | 1.68 | 1.67 ± 0.42 | 1000 |
| 4/41 | 30.90 | 30.7 ± 7.75 | 5.35 | 5.31 ± 1.34 | 1000 |
| 4/5 | 12.32 | 12.24 ± 3.09 | 69.15 | 68.7 ± 17.34 | 1000 |
Residual error was added to simulated data using a constant coefficient of variation structure, where the observed metabolic rate (Vobs)for the ithindividual at the jth substrate incubation concentration. The residual errors (εi) were normally distributed with a mean of 0 and a standard deviation (σ) set at the following values: 0, 0.05, 0.10, or 0.2; corresponding to a coefficient of variation (CV%) of 0%, 5%, 10%, and 20%.
# Error Design Datatable -------
error_design.data <- rbind(population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 0) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "0% CV"),
population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 5) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "5% CV"),
population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 10) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "10% CV"),
population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 20) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "20% CV")) %>%
filter(S %in% c(0.1, 0.5, 1, 2, 5, 10, 15, 25, 40, 55, 70, 85, 100)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Rich Design (9 pts)")
# Assigning Condition as Ordered Factor -----------
error_design.data$Condition <- factor(error_design.data$Condition,
levels = c("0% CV","5% CV","10% CV","20% CV"))
# Simulation Grid ---------
error.nls <- nls(formula = V ~ (Vmax*S)/(Km + S), data = error_design.data,
start = list(Vmax = 350, Km = 2))
error.grid <- tibble(S = seq(0,100,0.1),
V = (coef(error.nls)[[1]]*S)/(coef(error.nls)[[2]]+S))# Experimental Error Scenarios Plot -----------
ggplot(error_design.data %>% filter(S %in% c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)),
aes(x = S, y = V))+
geom_line(data = error.grid, size = 1, alpha = 0.75)+
geom_point(aes(y = V_adj, group = ID2, color = ID2, shape = ID2), size = 2.5)+
facet_wrap(~Condition, ncol = 2, scales = "free")+
theme_bw(base_size = 14)+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.7)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.9,0.11), legend.box = "horizontal",
legend.background = element_rect(color = "grey"))+
labs(title = "Experimental Error Scenarios")Illustration of experimental error scenarios for 3 subjects with the same genotype. Black line represents the non-linear least-squares Michaelis-Menten fit for the genotype group.
A custom function, update_data(), was created to
generate experimental data in-silico according to the specified
experimental design and conditions. Two in-silico experimental designs
(rich design, and sparse design) for two experimental conditions (UM/EM
and IM/PM conditions) were incorporated into the function. Rich design
experiments were conducted using 9 overlapping atomoxetine (ATX)
incubation concentrations per subject, while sparse design experiments
used 3-4 staggered incubation concentrations per subject , with both
designs covering similar concentration ranges according to their
experimental condition (UM/EM range: 0.10 - 100 uM, IM/PM range: 1-2000
uM). Additionally, each experimental design can subjected to variable
degrees of experimental error using the CV% input. The
inputs for the function the update_data() function include
the following:
n_indiv - number of subjects per diplotype group (default = 10).
CV%_D - Average standard deviation for each diplotype group as a percentage of the true estimate for Vmax and Km (between subject variability for each diplotype group) (default = 25).
CV% - Average coefficient of variation across all points for simulated Michaelis-Menten kinetics for each individual (residual error).
start - Starting concentration from Michaelis-Menten simulation (default = 0).
end - Ending concentration for Michaelis-Menten kinetics (default = 100 uM).
by - Simulated concentration increment (default = 0.5).
Pop_freq - TRUE or FALSE argument indicating whether to use true population frequency (default is FALSE).
Seed - Random number generator seed (for reproducibility, default = 23457).
type - Experimental design to be simulated (options: rich sampling (9pt) = “rich”, sparse sampling (3pt) = “sparse3pt”, sparse sampling (4pt) = “sparse4pt”, rich sampling for poor metabolizers (16 pt) = “PM”, sparse sampling for poor metabolizer (3pt) = “PM_3pt”, sparse sampling for poor metabolizer (4pt) = “PM_4pt”.
update_data <-function(`CV%`, type, n_indiv = 10, `CV%.D` = 25, by = 0.1,
start = 0, end = 100, Seed = 23457, Pop.freq = FALSE){
# Sampling Information --------------
# Full Range Set
full_set <- c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)
PM_range <- c(1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, 2000)
# Strategic Sampling Sets ----
sparse3pt_set1 <- c(0.1, 2.5, 50)
sparse3pt_set2 <- c(1, 10, 100)
sparse4pt_set1 <- c(0.1, 2.5, 25, 50)
sparse4pt_set2 <- c(1, 10, 50, 100)
PM_3pt_set1 <- c(10,100,1000)
PM_3pt_set2 <- c(25,250,2000)
PM_4pt_set1 <- c(1,10,100,1000)
PM_4pt_set2 <- c(5,25,250,2000)
PM_range_set1 <- c(1,10,25,50,100,400,1000,2000)
PM_range_set2 <- c(5,15,30,60,90,200,800,1600)
# Rich Sampling Simulation ----------------
Rich <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
# Sparse Sampling Simulation (3 point) ------------
Sparse3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse3pt_set1, S %in% sparse3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
# Sparse Sampling Simulation (4 point) --------------
Sparse4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse4pt_set1, S %in% sparse4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
PM <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
PM_3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_3pt_set1, S %in% PM_3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
PM_4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_4pt_set1, S %in% PM_4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
# Output Selection ----------
switch(type,
"rich" = Rich,
"sparse3pt" = Sparse3pt,
"sparse4pt" = Sparse4pt,
"PM" = PM,
"PM_3pt" = PM_3pt,
"PM_4pt" = PM_4pt)
}Strategic sampling was accomplished by evenly dividing individuals falling under the same CYP2D6 genotype into two groups (Group 1 and Group 2). Each group was designated a pre-defined set of incubation concentrations to be used for the in-silico experiment with the total range of incubation concentrations dependent on the type of CYP2D6 metabolizer (as described in the Strategic Sampling Approach section above). Note: strategic sampling was not used for rich sample designs, meaning all subject were subjected to the same set of incubations concentrations.
Extensive and Ultra-Rapid Metabolizers
Incubation Concentrations - Rich (9pt) - Set: 0.1, 0.5, 1, 2.5, 5, 10, 25, 50, and 100 uM
Incubation Concentrations - Sparse (3pt) - Set 1: 0.1, 2.5, and 50 uM - Set 2: 1, 10, and 100 uM
Incubation Concentrations - Sparse (4pt) - Set 1: 0.1, 2.5, 25, and 50 uM - Set 2: 1, 10, 50, and 100 uM
Intermediate and Poor Metabolizers
Incubation Concentrations - Rich (16pt) - Set: 1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, and 2000 uM
Incubation Concentrations - Sparse (3pt) - Set 1: 10,100, and 1000 uM - Set 2: 25, 250, and 2000 uM
Incubation Concentrations - Sparse (4pt) - Set 1: 1,10,100, and 1000 uM - Set 2: 5,25,250, and 2000 uM
# Generating Representative Data sets =========
## Representative data comes from 2 wild-type subjects.
rich_design.data <- update_data(`CV%` = 0, type = "rich") %>%
filter(Diplotype == "1/1", ID %in% c(1,2)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Rich Design (9 pts)")
sparse_3pt_design.data <- update_data(`CV%` = 0, type = "sparse3pt")%>%
filter(Diplotype == "1/1", ID %in% c(1,2)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (3 pts)")
sparse_4pt_design.data <- update_data(`CV%` = 0, type = "sparse4pt")%>%
filter(Diplotype == "1/1", ID %in% c(1,2)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (4 pts)")
PM_design.data <- update_data(`CV%` = 0, type = "PM") %>%
filter(Diplotype == "4/5", ID %in% c(81,82)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Rich Design (16 pts)")
PM3pt_design.data <- update_data(`CV%` = 0, type = "PM_3pt") %>%
filter(Diplotype == "4/5", ID %in% c(81,82)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (3 pts)")
PM4pt_design.data <- update_data(`CV%` = 0, type = "PM_4pt") %>%
filter(Diplotype == "4/5", ID %in% c(81,82)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (4 pts)")
# Solution Grid
design.nls <- nls(formula = V ~ (Vmax*S)/(Km + S), data = rich_design.data,
start = list(Vmax = 350, Km = 2))
grid <- tibble(S = seq(0,100,0.1),
V = (coef(design.nls)[[1]]*S)/(coef(design.nls)[[2]]+S))
design_PM.nls <- nls(formula = V ~ (Vmax*S)/(Km + S), data = PM_design.data,
start = list(Vmax = 20, Km = 50))
grid_pm <- tibble(S = seq(0,2000,0.1),
V = (coef(design_PM.nls)[[1]]*S)/(coef(design_PM.nls)[[2]]+S))# UM/EM STRATEGIC SAMPLING DESIGN FIGURE ================
## Representative plots contain data for 2 wild-type subjects.
# Rich Design--------
D1<- ggplot(rich_design.data, aes(x = S, y = V))+
geom_line(data = grid, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Conventional")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# Sparse Design (3pt)----------
D2<- ggplot(sparse_3pt_design.data, aes(x = S, y = V))+
geom_line(data = grid, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Strategic Scenario (2)")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# Sparse Design (4pt)------------
D3 <- ggplot(sparse_4pt_design.data, aes(x = S, y = V))+
geom_line(data = grid, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Strategic Scenario (1)")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# IM/PM STRATEGIC SAMPLING DESIGN FIGURE ================
## Representative plots contain data for 2 wild-type subjects.
# PM Sparse Design (4pt)---------
D4 <- ggplot(PM4pt_design.data, aes(x = S, y = V))+
geom_line(data = grid_pm, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# Extensize metabolizer model---------
D5 <- D3 + labs(title = "") #4pt
D6 <- D2 + labs(title = "") #3pt
D7 <- D1 + labs(title = "Extensive Metabolizer") #9pt
# PM Sparse Design (3pt)---------
D8 <- ggplot(PM3pt_design.data, aes(x = S, y = V))+
geom_line(data = grid_pm, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# PM Rich Design (16pt)---------
D9 <- ggplot(PM_design.data, aes(x = S, y = V))+
geom_line(data = grid_pm, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Poor Metabolizers")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)## Extensive Metabolizers
ggpubr::ggarrange(D7,D5,D6,ncol = 2, nrow = 2, labels = "AUTO")## Poor Metabolizers
ggpubr::ggarrange(D9,D4,D8, ncol = 2, nrow = 2, labels = "AUTO")# Population Datasets (CV = 0%)
rich.data <- update_data(`CV%` = 0, type = "rich") %>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt.data <- update_data(`CV%` = 0, type = "sparse3pt") %>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt.data <- update_data(`CV%` = 0, type = "sparse4pt") %>%
filter(!Diplotype %in% c("4/41", "4/5"))A grouped data frame was created for each experiment where reaction
rate (V) is predicted by substrate concentration (S) according to
subject (ID). Individual fits (un-weighted) and parameter estimates for
the data sets were obtained using non-linear least squares function
nlsList(). Values obtained from the nlsList()
were used as initial estimates for the non-linear mixed effect
model.
# Grouped Data frame -------------
rich_data.grp <- groupedData(V~S|ID, rich.data)
sparse3pt.grp <- groupedData(V~S|ID, sparse3pt.data)
sparse4pt.grp <- groupedData(V~S|ID, sparse4pt.data)
# Fit Model ((Non-Linear Least Squares) ------------
rich_fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_data.grp)
sparse3pt_fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt.grp)
sparse4pt_fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt.grp)Non-linear mixed effect modeling was conducted using the
nlme() package version 3.1 (Bates et al., 2000). The
nlsList() fits were used as initial estimates for the NLME
model. Individual level random effects were applied to both Vmax and Km.
A diagonal variance-covariance structure was assigned to each model
using the pdDiag() function.
# Fit Model ((Non-Linear Mixed-Effect) ------------
rich.nlme <- nlme(rich_fit.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt.nlme <- nlme(sparse3pt_fit.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt.nlme <- nlme(sparse4pt_fit.nls, random = pdDiag(Vmax + Km ~ 1))
# Model Fixed Effect comparisons
fixef.table <- rbind(fixef(rich.nlme),
fixef(sparse3pt.nlme),
fixef(sparse4pt.nlme)) %>%
as_tibble() %>%
mutate(Model = c("Rich","Sparse (3pt)","Sparse (4pt)")) %>%
relocate(Model)
fixef.table %>%
mutate(across(where(is.numeric), round, 2)) %>%
kbl( table.attr = "style='width:60%;'",
caption = "Base Model: Mixed Effect Michaelis-Menten") %>%
kable_classic("striped")| Model | Vmax | Km |
|---|---|---|
| Rich | 459.63 | 1.95 |
| Sparse (3pt) | 459.63 | 1.95 |
| Sparse (4pt) | 459.63 | 1.95 |
CYP2D6 genotype was incorporated as a categorical covariate in the model by adding it as a fixed effect that impacts both Vmax and Km. The extracted population level fixed effects for Vmax and Km from the base model (described in the previous section) were used as initial estimates for the covariate model.
#Models with genotype as a Covariate -----
rich.fxf <- fixef.table[1,2:3]
sparse3pt.fxf <- fixef.table[2,2:3]
sparse4pt.fxf <- fixef.table[3,2:3]
rich_covar.nlme <- update(rich.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich.fxf[[1]], rep(0,6),
rich.fxf[[2]], rep(0,6)))
sparse3pt_covar.nlme <- update(sparse3pt.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt.fxf[[1]], rep(0,6),
sparse3pt.fxf[[2]], rep(0,6)))
sparse4pt_covar.nlme <- update(sparse4pt.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt.fxf[[1]], rep(0,6),
sparse4pt.fxf[[2]], rep(0,6)))The impact of adding CYP2D6 genotype as a (covariate on both Vmax and
Km) and the whether adding covariates to the model improved fit compared
to baseline was assessed using analysis of variance
anova().
anova(rich.nlme, rich_covar.nlme)
anova(sparse3pt.nlme, sparse3pt_covar.nlme)
anova(sparse4pt.nlme, sparse4pt_covar.nlme)anova(rich_covar.nlme) %>%
kbl(caption = "Rich Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")
anova(sparse3pt_covar.nlme)%>%
kbl(caption = "Sparse 3pt Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")
anova(sparse4pt_covar.nlme)%>%
kbl(caption = "Sparse 4pt Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")# Covariate Analysis Summary Table ------------
covaref.table <- rbind(fixef(rich_covar.nlme),
fixef(sparse3pt_covar.nlme),
fixef(sparse4pt_covar.nlme)) %>%
as_tibble() %>%
mutate(`Covariate Model` = c("Rich","Sparse (3pt)","Sparse (4pt)")) %>%
pivot_longer(`Vmax.(Intercept)`:`Km.Diplotype2/4`,
names_to = "Parameter",
values_to = "Value") %>%
pivot_wider(names_from = `Covariate Model`, values_from = Value) %>%
separate(Parameter, sep = "\\.", into = c("Parameter","Diplotype"), remove = T) %>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (1/1)"),
across(where(is.numeric), round, 2))# Reference Estimates -----------
covaref.table %>%
filter(Diplotype == "Reference (1/1)")%>%
kbl(caption = "Reference Population Estimates",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")| Parameter | Diplotype | Rich | Sparse (3pt) | Sparse (4pt) |
|---|---|---|---|---|
| Vmax | Reference (1/1) | 367.01 | 367.01 | 367.01 |
| Km | Reference (1/1) | 1.50 | 1.50 | 1.50 |
# CYP2D6 Genotype Effects --------
covaref.table%>%
kbl(caption = "Michaelis Menten Mixed-Effect Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped") %>%
pack_rows("Vmax Estimate (Population)", 1,1) %>%
pack_rows("Covariate (Vmax Eta)", 2, 7)%>%
pack_rows("Km Estimate (Population)", 8,8) %>%
pack_rows("Covariate (Km Eta)", 9, 14)| Parameter | Diplotype | Rich | Sparse (3pt) | Sparse (4pt) |
|---|---|---|---|---|
| Vmax Estimate (Population) | ||||
| Vmax | Reference (1/1) | 367.01 | 367.01 | 367.01 |
| Covariate (Vmax Eta) | ||||
| Vmax | Diplotype1/2 | 119.92 | 119.92 | 119.92 |
| Vmax | Diplotype1/2x2 | 1092.74 | 1092.74 | 1092.74 |
| Vmax | Diplotype1/4 | -83.34 | -83.34 | -83.34 |
| Vmax | Diplotype2/2 | -106.12 | -106.11 | -106.12 |
| Vmax | Diplotype2/3 | -107.15 | -107.15 | -107.15 |
| Vmax | Diplotype2/4 | -267.73 | -267.72 | -267.73 |
| Km Estimate (Population) | ||||
| Km | Reference (1/1) | 1.50 | 1.50 | 1.50 |
| Covariate (Km Eta) | ||||
| Km | Diplotype1/2 | -0.06 | -0.06 | -0.06 |
| Km | Diplotype1/2x2 | 0.39 | 0.39 | 0.39 |
| Km | Diplotype1/4 | -0.55 | -0.55 | -0.55 |
| Km | Diplotype2/2 | 3.25 | 3.25 | 3.25 |
| Km | Diplotype2/3 | -0.13 | -0.13 | -0.13 |
| Km | Diplotype2/4 | 0.23 | 0.23 | 0.23 |
# Function to Extract Confidence Interval Info ----------
extract_CI <- function(int.data, name){
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(across(where(is.numeric), round, 2),
`CI (95%)` = paste0("(",lower,", ",upper,")")) %>%
select(est., `CI (95%)`)
colnames(temp) <- c(paste(name, " Estimate"), paste(name, " CI (95%)"))
temp
}
# Confidence Interval by Sampling Scenario -------
rich_covar.intervals <- intervals(rich_covar.nlme, which = "fixed")
sparse3pt_covar.intervals <- intervals(sparse3pt_covar.nlme, which = "fixed")
sparse4pt_covar.intervals <- intervals(sparse4pt_covar.nlme, which = "fixed")
# Summary Table of Confidence Intervals by Diplotype and Scenario -----------
covar_ci.table <- cbind(extract_CI(rich_covar.intervals, name = "Rich"),
extract_CI(sparse3pt_covar.intervals, name = "Sparse 3pt"),
extract_CI(sparse4pt_covar.intervals, name = "Sparse 4pt")) %>%
mutate(Parameter = rownames(rich_covar.intervals$fixed)) %>%
relocate(Parameter) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (1/1)"))
colnames(covar_ci.table)<- c("Parameter", "Diplotype",
"Estimate","CI (95%)",
"Estimate", "CI (95%)",
"Estimate", "CI (95%)")
covar_ci.table%>%
kbl(caption = "Michaelis Menten Mixed-Effect Model (w/Covariates)") %>%
kable_classic(full_width = T, "striped")%>%
pack_rows("Vmax Estimate (Population)", 1,1) %>%
pack_rows("Covariate (Vmax Eta)", 2, 7)%>%
pack_rows("Km Estimate (Population)", 8,8) %>%
pack_rows("Covariate (Km Eta)", 9, 14) %>%
add_header_above(c(" ", " ",
"Rich (9pt)" = 2,
"Sparse (3pt)" = 2,
"Sparse (4pt)" = 2))|
Rich (9pt)
|
Sparse (3pt)
|
Sparse (4pt)
|
|||||
|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Estimate | CI (95%) | Estimate | CI (95%) | Estimate | CI (95%) |
| Vmax Estimate (Population) | |||||||
| Vmax | Reference (1/1) | 367.01 | (319.69, 414.33) | 367.01 | (319.34, 414.67) | 367.01 | (319.51, 414.51) |
| Covariate (Vmax Eta) | |||||||
| Vmax | Diplotype1/2 | 119.92 | (53.01, 186.84) | 119.92 | (52.51, 187.33) | 119.92 | (52.74, 187.1) |
| Vmax | Diplotype1/2x2 | 1092.74 | (1025.83, 1159.66) | 1092.74 | (1025.33, 1160.16) | 1092.74 | (1025.56, 1159.93) |
| Vmax | Diplotype1/4 | -83.34 | (-150.26, -16.42) | -83.34 | (-150.75, -15.93) | -83.34 | (-150.52, -16.16) |
| Vmax | Diplotype2/2 | -106.12 | (-173.03, -39.2) | -106.11 | (-173.52, -38.7) | -106.12 | (-173.3, -38.94) |
| Vmax | Diplotype2/3 | -107.15 | (-174.07, -40.24) | -107.15 | (-174.56, -39.74) | -107.15 | (-174.33, -39.97) |
| Vmax | Diplotype2/4 | -267.73 | (-334.64, -200.81) | -267.72 | (-335.13, -200.31) | -267.73 | (-334.91, -200.54) |
| Km Estimate (Population) | |||||||
| Km | Reference (1/1) | 1.50 | (1.33, 1.67) | 1.50 | (1.33, 1.67) | 1.50 | (1.33, 1.67) |
| Covariate (Km Eta) | |||||||
| Km | Diplotype1/2 | -0.06 | (-0.31, 0.18) | -0.06 | (-0.31, 0.18) | -0.06 | (-0.31, 0.18) |
| Km | Diplotype1/2x2 | 0.39 | (0.14, 0.63) | 0.39 | (0.14, 0.63) | 0.39 | (0.14, 0.63) |
| Km | Diplotype1/4 | -0.55 | (-0.79, -0.3) | -0.55 | (-0.79, -0.3) | -0.55 | (-0.79, -0.3) |
| Km | Diplotype2/2 | 3.25 | (3.01, 3.49) | 3.25 | (3.01, 3.5) | 3.25 | (3.01, 3.5) |
| Km | Diplotype2/3 | -0.13 | (-0.38, 0.11) | -0.13 | (-0.38, 0.11) | -0.13 | (-0.38, 0.11) |
| Km | Diplotype2/4 | 0.23 | (-0.01, 0.48) | 0.23 | (-0.01, 0.48) | 0.23 | (-0.01, 0.48) |
Data sets which include increasing degrees of residual error (5, 10,
and 20% CV) were generated using the update_data() function
described earlier (see In-Silico Experiments section for full
description). Additonally, the same methodology for model analysis used
in the the previous section (Population Michaelis-Menten Modeling
(UM/EM)) was applied here to evaluate the impact of increasing
experimental error on parameter estimates.
# Population Data sets (CV = 5%) ----------------
rich_CV5.data <- update_data(`CV%` = 5, type = "rich")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt_CV5.data <- update_data(`CV%` = 5, type = "sparse3pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt_CV5.data <- update_data(`CV%` = 5, type = "sparse4pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
# Population Data sets (CV = 10%) ----------------
rich_CV10.data <- update_data(`CV%` = 10, type = "rich")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt_CV10.data <- update_data(`CV%` = 10, type = "sparse3pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt_CV10.data <- update_data(`CV%` = 10, type = "sparse4pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
# Population Data sets (CV = 20%) ----------------
rich_CV20.data <- update_data(`CV%` = 20, type = "rich")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt_CV20.data <- update_data(`CV%` = 20, type = "sparse3pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt_CV20.data <- update_data(`CV%` = 20, type = "sparse4pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))# Grouped Data frame =============================
rich_CV5_data.grp <- groupedData(V~S|ID, rich_CV5.data)
sparse3pt_CV5.grp <- groupedData(V~S|ID, sparse3pt_CV5.data)
sparse4pt_CV5.grp <- groupedData(V~S|ID, sparse4pt_CV5.data)
rich_CV10_data.grp <- groupedData(V~S|ID, rich_CV10.data)
sparse3pt_CV10.grp <- groupedData(V~S|ID, sparse3pt_CV10.data)
sparse4pt_CV10.grp <- groupedData(V~S|ID, sparse4pt_CV10.data)
rich_CV20_data.grp <- groupedData(V~S|ID, rich_CV20.data)
sparse3pt_CV20.grp <- groupedData(V~S|ID, sparse3pt_CV20.data)
sparse4pt_CV20.grp <- groupedData(V~S|ID, sparse4pt_CV20.data)
# Fit Model ((Non-Linear Least Squares) =====================
rich_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_CV5_data.grp)
sparse3pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt_CV5.grp)
sparse4pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt_CV5.grp)
rich_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_CV10_data.grp)
sparse3pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt_CV10.grp)
sparse4pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt_CV10.grp)
rich_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_CV20_data.grp)
sparse3pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt_CV20.grp)
sparse4pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt_CV20.grp)# Fit Model ((Non-Linear Mixed-Effect) ===================
rich_CV5.nlme <- nlme(rich_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt_CV5.nlme <- nlme(sparse3pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt_CV5.nlme <- nlme(sparse4pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
rich_CV10.nlme <- nlme(rich_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt_CV10.nlme <- nlme(sparse3pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt_CV10.nlme <- nlme(sparse4pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
rich_CV20.nlme <- nlme(rich_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt_CV20.nlme <- nlme(sparse3pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt_CV20.nlme <- nlme(sparse4pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))# Model Fixed Effect comparisons----------
fixef_CV.table <- rbind(
fixef(rich.nlme),
fixef(sparse3pt.nlme),
fixef(sparse4pt.nlme),
fixef(rich_CV5.nlme),
fixef(sparse3pt_CV5.nlme),
fixef(sparse4pt_CV5.nlme),
fixef(rich_CV10.nlme),
fixef(sparse3pt_CV10.nlme),
fixef(sparse4pt_CV10.nlme),
fixef(rich_CV20.nlme),
fixef(sparse3pt_CV20.nlme),
fixef(sparse4pt_CV20.nlme)) %>%
as_tibble() %>%
mutate(Model = c("Rich (0% CV)","Sparse (0% CV, 3pt)","Sparse (0% CV, 4pt)",
"Rich (5% CV)","Sparse (5% CV, 3pt)","Sparse (5% CV, 4pt)",
"Rich (10% CV)","Sparse (10% CV, 3pt)","Sparse (10% CV, 4pt)",
"Rich (20% CV)","Sparse (20% CV, 3pt)","Sparse (20% CV, 4pt)"),
across(where(is.numeric),round,3)) %>%
relocate(Model)
# Fixed Effect Table across conditions -------------------
fixef_CV.table %>%
kbl(caption = "Population Estimates (w/ Residual Error)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")%>%
pack_rows(" ", 1,3) %>%
pack_rows(" ", 4, 6)%>%
pack_rows(" ", 7, 9)%>%
pack_rows(" ", 10, 12)| Model | Vmax | Km |
|---|---|---|
| Rich (0% CV) | 459.628 | 1.948 |
| Sparse (0% CV, 3pt) | 459.628 | 1.948 |
| Sparse (0% CV, 4pt) | 459.628 | 1.948 |
| Rich (5% CV) | 458.428 | 1.689 |
| Sparse (5% CV, 3pt) | 460.684 | 1.726 |
| Sparse (5% CV, 4pt) | 459.051 | 1.692 |
| Rich (10% CV) | 461.150 | 1.723 |
| Sparse (10% CV, 3pt) | 465.169 | 1.888 |
| Sparse (10% CV, 4pt) | 461.482 | 1.804 |
| Rich (20% CV) | 464.268 | 1.765 |
| Sparse (20% CV, 3pt) | 468.840 | 1.924 |
| Sparse (20% CV, 4pt) | 464.632 | 1.924 |
# Function to extract confidence intervals from model ---------
extract_CI2 <- function(int.data, name){
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(across(where(is.numeric), round, 2),
`CI (95%)` = paste0("(",lower,", ",upper,")")) %>%
select(est., `CI (95%)`)
colnames(temp) <- c("Estimate","CI (95%)")
temp$Model <- name
temp$Parameter <- c("Vmax","Km")
temp
}
# Confidence Interval Data ---------
rich_CV5.intervals <- intervals(rich_CV5.nlme, which = "fixed")
sparse3pt_CV5.intervals <- intervals(sparse3pt_CV5.nlme, which = "fixed")
sparse4pt_CV5.intervals <- intervals(sparse4pt_CV5.nlme, which = "fixed")
rich_CV10.intervals <- intervals(rich_CV10.nlme, which = "fixed")
sparse3pt_CV10.intervals <- intervals(sparse3pt_CV10.nlme, which = "fixed")
sparse4pt_CV10.intervals <- intervals(sparse4pt_CV10.nlme, which = "fixed")
rich_CV20.intervals <- intervals(rich_CV20.nlme, which = "fixed")
sparse3pt_CV20.intervals <- intervals(sparse3pt_CV20.nlme, which = "fixed")
sparse4pt_CV20.intervals <- intervals(sparse4pt_CV20.nlme, which = "fixed")
# Confidence Interval Tables ----------
CV5_ci.table <- rbind(extract_CI2(rich_CV5.intervals, name = "Rich (5% CV)"),
extract_CI2(sparse3pt_CV5.intervals, name = "Sparse 3pt (5% CV)"),
extract_CI2(sparse4pt_CV5.intervals, name = "Sparse 4pt (5% CV)")) %>%
relocate(Model, Parameter) %>%
arrange(desc(Parameter))
CV10_ci.table <- rbind(extract_CI2(rich_CV10.intervals, name = "Rich (10% CV)"),
extract_CI2(sparse3pt_CV10.intervals, name = "Sparse 3pt (10% CV)"),
extract_CI2(sparse4pt_CV10.intervals, name = "Sparse 4pt (10% CV)")) %>%
relocate(Model, Parameter) %>%
arrange(desc(Parameter))
CV20_ci.table <- rbind(extract_CI2(rich_CV20.intervals, name = "Rich (20% CV)"),
extract_CI2(sparse3pt_CV20.intervals, name = "Sparse 3pt (20% CV)"),
extract_CI2(sparse4pt_CV20.intervals, name = "Sparse 4pt (20% CV)")) %>%
relocate(Model, Parameter) %>%
arrange(desc(Parameter))
# Summary Table of Estimates with CIs ----------
rbind(CV5_ci.table, CV10_ci.table, CV20_ci.table) %>%
arrange(desc(Parameter))%>%
kbl(caption = "Population Estimates (w/ Residual Error & 95% CI)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")%>%
pack_rows("Vmax Estimates", 1,9) %>%
pack_rows("", 1, 3) %>%
pack_rows("", 4, 6) %>%
pack_rows("", 7, 9) %>%
pack_rows("Km Estimates", 10, 18) %>%
pack_rows(" ", 10, 12) %>%
pack_rows(" ", 13, 15) %>%
pack_rows(" ", 16, 18)| Model | Parameter | Estimate | CI (95%) |
|---|---|---|---|
| Vmax Estimates | |||
| Rich (5% CV) | Vmax | 458.43 | (356.53, 560.33) |
| Sparse 3pt (5% CV) | Vmax | 460.68 | (358.37, 563) |
| Sparse 4pt (5% CV) | Vmax | 459.05 | (357.03, 561.07) |
| Rich (10% CV) | Vmax | 461.15 | (358.84, 563.46) |
| Sparse 3pt (10% CV) | Vmax | 465.17 | (362.37, 567.97) |
| Sparse 4pt (10% CV) | Vmax | 461.48 | (358.88, 564.08) |
| Rich (20% CV) | Vmax | 464.27 | (360.44, 568.1) |
| Sparse 3pt (20% CV) | Vmax | 468.84 | (363.07, 574.61) |
| Sparse 4pt (20% CV) | Vmax | 464.63 | (360.16, 569.1) |
| Km Estimates | |||
| Rich (5% CV) | Km | 1.69 | (1.52, 1.86) |
| Sparse 3pt (5% CV) | Km | 1.73 | (1.59, 1.86) |
| Sparse 4pt (5% CV) | Km | 1.69 | (1.55, 1.83) |
| Rich (10% CV) | Km | 1.72 | (1.59, 1.85) |
| Sparse 3pt (10% CV) | Km | 1.89 | (1.74, 2.03) |
| Sparse 4pt (10% CV) | Km | 1.80 | (1.62, 1.99) |
| Rich (20% CV) | Km | 1.77 | (1.61, 1.93) |
| Sparse 3pt (20% CV) | Km | 1.92 | (1.66, 2.19) |
| Sparse 4pt (20% CV) | Km | 1.92 | (1.66, 2.19) |
# Models with Diplotype as a Covariate -----
rich_CV5.fxf <- fixef_CV.table[1,2:3]
sparse3pt_CV5.fxf <- fixef_CV.table[2,2:3]
sparse4pt_CV5.fxf <- fixef_CV.table[3,2:3]
rich_CV10.fxf <- fixef_CV.table[4,2:3]
sparse3pt_CV10.fxf <- fixef_CV.table[5,2:3]
sparse4pt_CV10.fxf <- fixef_CV.table[6,2:3]
rich_CV20.fxf <- fixef_CV.table[7,2:3]
sparse3pt_CV20.fxf <- fixef_CV.table[8,2:3]
sparse4pt_CV20.fxf <- fixef_CV.table[9,2:3]
# NLME Covariate Model w/5% CV ==================
rich_CV5_covar.nlme <- update(rich_CV5.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich_CV5.fxf[[1]], rep(0,6),
rich_CV5.fxf[[2]], rep(0,6)))
sparse3pt_CV5_covar.nlme <- update(sparse3pt_CV5.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt_CV5.fxf[[1]], rep(0,6),
sparse3pt_CV5.fxf[[2]], rep(0,6)))
sparse4pt_CV5_covar.nlme <- update(sparse4pt_CV5.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt_CV5.fxf[[1]], rep(0,6),
sparse4pt_CV5.fxf[[2]], rep(0,6)))
# NLME Covariate Model w/10% CV ==================
rich_CV10_covar.nlme <- update(rich_CV10.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich_CV10.fxf[[1]], rep(0,6),
rich_CV10.fxf[[2]], rep(0,6)))
sparse3pt_CV10_covar.nlme <- update(sparse3pt_CV10.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt_CV10.fxf[[1]], rep(0,6),
sparse3pt_CV10.fxf[[2]], rep(0,6)))
sparse4pt_CV10_covar.nlme <- update(sparse4pt_CV10.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt_CV10.fxf[[1]], rep(0,6),
sparse4pt_CV10.fxf[[2]], rep(0,6)))
# NLME Covariate Model w/20% CV ==================
rich_CV20_covar.nlme <- update(rich_CV20.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich_CV20.fxf[[1]], rep(0,6),
rich_CV20.fxf[[2]], rep(0,6)))
sparse3pt_CV20_covar.nlme <- update(sparse3pt_CV20.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt_CV20.fxf[[1]], rep(0,6),
sparse3pt_CV20.fxf[[2]], rep(0,6)))
sparse4pt_CV20_covar.nlme <- update(sparse4pt_CV20.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt_CV20.fxf[[1]], rep(0,6),
sparse4pt_CV20.fxf[[2]], rep(0,6)))# Summary of Covariate Models (w/ variability) -------
anova(rich_CV5_covar.nlme) # RICH SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 547 2321.6898 <.0001
## Vmax.Diplotype 6 547 370.0573 <.0001
## Km.(Intercept) 1 547 2549.1590 <.0001
## Km.Diplotype 6 547 30.1046 <.0001
anova(sparse3pt_CV5_covar.nlme) # SPARSE 3PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 127 2590.3045 <.0001
## Vmax.Diplotype 6 127 366.3843 <.0001
## Km.(Intercept) 1 127 2536.6700 <.0001
## Km.Diplotype 6 127 19.4802 <.0001
anova(sparse4pt_CV5_covar.nlme) # SPARSE 4PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 197 2515.2320 <.0001
## Vmax.Diplotype 6 197 357.4567 <.0001
## Km.(Intercept) 1 197 2596.3871 <.0001
## Km.Diplotype 6 197 19.9836 <.0001
anova(rich_CV10_covar.nlme) # RICH SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 547 2171.6492 <.0001
## Vmax.Diplotype 6 547 310.4159 <.0001
## Km.(Intercept) 1 547 1763.7233 <.0001
## Km.Diplotype 6 547 11.3821 <.0001
anova(sparse3pt_CV10_covar.nlme) # SPARSE 3PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 127 2382.5424 <.0001
## Vmax.Diplotype 6 127 340.1663 <.0001
## Km.(Intercept) 1 127 756.2901 <.0001
## Km.Diplotype 6 127 6.9697 <.0001
anova(sparse4pt_CV10_covar.nlme) # SPARSE 4PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 197 2366.5954 <.0001
## Vmax.Diplotype 6 197 339.9690 <.0001
## Km.(Intercept) 1 197 741.8162 <.0001
## Km.Diplotype 6 197 6.7521 <.0001
anova(rich_CV20_covar.nlme) # RICH SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 547 1768.6179 <.0001
## Vmax.Diplotype 6 547 258.1663 <.0001
## Km.(Intercept) 1 547 457.3300 <.0001
## Km.Diplotype 6 547 2.9774 0.0072
anova(sparse3pt_CV20_covar.nlme) # SPARSE 3PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 127 1490.6401 <.0001
## Vmax.Diplotype 6 127 217.4530 <.0001
## Km.(Intercept) 1 127 193.1641 <.0001
## Km.Diplotype 6 127 2.7206 0.0161
anova(sparse4pt_CV20_covar.nlme) # SPARSE 4PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 197 1716.8079 <.0001
## Vmax.Diplotype 6 197 252.6265 <.0001
## Km.(Intercept) 1 197 189.6812 <.0001
## Km.Diplotype 6 197 2.5226 0.0225
Heterogeneity of variance was examined qualitatively by plotting
standardized residules (Pearson). Heteroscedacity was corrected by
modeleling residual variance as a power function of mean fitted value
using the varPower() function.
plot(rich_covar.nlme, main = "Rich (9pt, 0% CV)")plot(rich_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 0% CV, weighted)")plot(sparse3pt_covar.nlme, main = "Sparse (3pt, 0% CV)")plot(sparse3pt_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 0% CV, weighted)")plot(sparse4pt_covar.nlme, main = "Sparse (4pt, 0% CV)")plot(sparse4pt_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 0% CV, weighted)")plot(rich_CV5_covar.nlme, main = "Rich (9pt, 5% CV)")plot(rich_CV5_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 5% CV, weighted)")plot(sparse3pt_CV5_covar.nlme, main = "Sparse (3pt, 5% CV)")plot(sparse3pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 5% CV, weighted)")plot(sparse4pt_CV5_covar.nlme, main = "Sparse (4pt, 5% CV)")plot(sparse4pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 5% CV, weighted)")plot(rich_CV10_covar.nlme, main = "Rich (9pt, 10% CV)")plot(rich_CV10_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 10% CV, weighted)")plot(sparse3pt_CV10_covar.nlme, main = "Sparse (3pt, 10% CV)")plot(sparse3pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 10% CV, weighted)")plot(sparse4pt_CV10_covar.nlme, main = "Sparse (4pt, 10% CV)")plot(sparse4pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 10% CV, weighted)")plot(rich_CV20_covar.nlme, main = "Rich (9pt, 20% CV)")plot(rich_CV20_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 20% CV, weighted)")plot(sparse3pt_CV20_covar.nlme, main = "Sparse (3pt, 20% CV)")plot(sparse3pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 20% CV, weighted)")plot(sparse4pt_CV20_covar.nlme, main = "Sparse (4pt, 20% CV)")plot(sparse4pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 20% CV, weighted)")# Adjusting Residual Error Model ----
rich_CV5_covar.nlme <- rich_CV5_covar.nlme %>% update(weights = varPower())
sparse3pt_CV5_covar.nlme <- sparse3pt_CV5_covar.nlme %>% update(weights = varPower())
sparse4pt_CV5_covar.nlme <- sparse4pt_CV5_covar.nlme %>% update(weights = varPower())
rich_CV10_covar.nlme <- rich_CV10_covar.nlme %>% update(weights = varPower())
sparse3pt_CV10_covar.nlme <- sparse3pt_CV10_covar.nlme %>% update(weights = varPower())
sparse4pt_CV10_covar.nlme <- sparse4pt_CV10_covar.nlme %>% update(weights = varPower())
rich_CV20_covar.nlme <- rich_CV20_covar.nlme %>% update(weights = varPower())
sparse3pt_CV20_covar.nlme <- sparse3pt_CV20_covar.nlme %>% update(weights = varPower())
sparse4pt_CV20_covar.nlme <- sparse4pt_CV20_covar.nlme %>% update(weights = varPower())# Creating covariate fixed effect table ------
covaref_CV.table <- rbind(fixef(rich_CV5_covar.nlme),
fixef(sparse3pt_CV5_covar.nlme),
fixef(sparse4pt_CV5_covar.nlme),
fixef(rich_CV10_covar.nlme),
fixef(sparse3pt_CV10_covar.nlme),
fixef(sparse4pt_CV10_covar.nlme),
fixef(rich_CV20_covar.nlme),
fixef(sparse3pt_CV20_covar.nlme),
fixef(sparse4pt_CV20_covar.nlme)) %>%
as_tibble() %>%
mutate(`Covariate Model` = c("Rich (5% CV)","Sparse (5% CV, 3pt)","Sparse (5% CV, 4pt)",
"Rich (10% CV)","Sparse (10% CV, 3pt)","Sparse (10% CV, 4pt)",
"Rich (20% CV)","Sparse (20% CV, 3pt)","Sparse (20% CV, 4pt)")) %>%
pivot_longer(`Vmax.(Intercept)`:`Km.Diplotype2/4`,
names_to = "Parameter",
values_to = "Value") %>%
pivot_wider(names_from = `Covariate Model`, values_from = Value) %>%
separate(Parameter, sep = "\\.", into = c("Parameter","Diplotype"), remove = T) %>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (1/1)"),
across(where(is.numeric), round, 2))
# Complete Covariate Effects Table --------------------
covar_Vmax.table <- cbind(covaref.table, covaref_CV.table %>%
select(-Parameter, -Diplotype)) %>% filter(Parameter == "Vmax")
covar_Km.table <- cbind(covaref.table, covaref_CV.table %>%
select(-Parameter, -Diplotype)) %>% filter(Parameter == "Km")
# Vmax Estimates -------------------------------
final_vmax.est <- covar_Vmax.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (1/1)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (1/1)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (1/1)" = "Diplotype1/1")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Km), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Vmax)/Vmax)*100, digits = 2),
Vmax = round(Vmax, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Vmax, Predicted, `Error (%)`) %>%
rename("Predicted (Vmax)" = "Predicted") %>%
arrange(Condition)
# Km Estimates -----------------------------------
final_km.est <- covar_Km.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (1/1)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (1/1)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (1/1)" = "Diplotype1/1")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Vmax), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Km)/Km)*100, digits = 2),
Km = round(Km, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Km, Predicted, `Error (%)`) %>%
rename("Predicted (Km)" = "Predicted") %>%
arrange(Condition)# Full Population Simulation rich at 0% CV ------------
# Rich Sampling (16 pt, PM) ---------------
PM.data <- update_data(n_indiv = 10, `CV%` = 0, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM.grp <- groupedData(V~S|ID, PM.data)
PM.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM.grp)
PM.nlme <- nlme(PM.nls, random = pdDiag(Vmax + Km ~ 1))
PM_covar.nlme <- update(PM.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM.nlme)[[1]], rep(0,1),
fixed.effects(PM.nlme)[[2]], rep(0,1)))
# Sparse Sampling (3pt, PM) -------------
PM_3pt.data <- update_data(n_indiv = 10, `CV%` = 0, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt.grp <- groupedData(V~S|ID, PM_3pt.data)
PM_3pt.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt.grp)
PM_3pt.nlme <- nlme(PM_3pt.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_covar.nlme <- update(PM_3pt.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt.nlme)[[2]], rep(0,1)))
# Sparse Sampling (4pt, PM) ---------------
PM_4pt.data <- update_data(n_indiv = 10, `CV%` = 0, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt.grp <- groupedData(V~S|ID, PM_4pt.data)
PM_4pt.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt.grp)
PM_4pt.nlme <- nlme(PM_4pt.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_covar.nlme <- update(PM_4pt.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt.nlme)[[2]], rep(0,1)))# 5% Residual Error ---------------------------
PM_CV5.data <- update_data(n_indiv = 10, `CV%` = 5, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_CV5_data.grp <- groupedData(V~S|ID, PM_CV5.data)
PM_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_CV5_data.grp)
PM_CV5.nlme <- nlme(PM_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
PM_CV5_covar.nlme <- update(PM_CV5.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_CV5.nlme)[[1]], rep(0,1),
fixed.effects(PM_CV5.nlme)[[2]], rep(0,1)))
PM_3pt_CV5.data <- update_data(n_indiv = 10, `CV%` = 5, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt_CV5_data.grp <- groupedData(V~S|ID, PM_3pt_CV5.data)
PM_3pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt_CV5_data.grp)
PM_3pt_CV5.nlme <- nlme(PM_3pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_CV5_covar.nlme <- update(PM_3pt_CV5.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt_CV5.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt_CV5.nlme)[[2]], rep(0,1)))
PM_4pt_CV5.data <- update_data(n_indiv = 10, `CV%` = 5, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt_CV5_data.grp <- groupedData(V~S|ID, PM_4pt_CV5.data)
PM_4pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt_CV5_data.grp)
PM_4pt_CV5.nlme <- nlme(PM_4pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_CV5_covar.nlme <- update(PM_4pt_CV5.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt_CV5.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt_CV5.nlme)[[2]], rep(0,1)))
# 10% Residual Error --------------------
PM_CV10.data <- update_data(n_indiv = 10, `CV%` = 10, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_CV10_data.grp <- groupedData(V~S|ID, PM_CV10.data)
PM_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_CV10_data.grp)
PM_CV10.nlme <- nlme(PM_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
PM_CV10_covar.nlme <- update(PM_CV10.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_CV10.nlme)[[1]], rep(0,1),
fixed.effects(PM_CV10.nlme)[[2]], rep(0,1)))
PM_3pt_CV10.data <- update_data(n_indiv = 10, `CV%` = 10, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt_CV10_data.grp <- groupedData(V~S|ID, PM_3pt_CV10.data)
PM_3pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt_CV10_data.grp)
PM_3pt_CV10.nlme <- nlme(PM_3pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_CV10_covar.nlme <- update(PM_3pt_CV10.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt_CV10.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt_CV10.nlme)[[2]], rep(0,1)))
PM_4pt_CV10.data <- update_data(n_indiv = 10, `CV%` = 10, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt_CV10_data.grp <- groupedData(V~S|ID, PM_4pt_CV10.data)
PM_4pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt_CV10_data.grp)
PM_4pt_CV10.nlme <- nlme(PM_4pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_CV10_covar.nlme <- update(PM_4pt_CV10.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt_CV10.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt_CV10.nlme)[[2]], rep(0,1)))
# 20% Residual Error ----------
PM_CV20.data <- update_data(n_indiv = 10, `CV%` = 20, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_CV20_data.grp <- groupedData(V~S|ID, PM_CV20.data)
PM_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_CV20_data.grp)
PM_CV20.nlme <- nlme(PM_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
PM_CV20_covar.nlme <- update(PM_CV20.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_CV20.nlme)[[1]], rep(0,1),
fixed.effects(PM_CV20.nlme)[[2]], rep(0,1)))
PM_3pt_CV20.data <- update_data(n_indiv = 10, `CV%` = 20, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt_CV20_data.grp <- groupedData(V~S|ID, PM_3pt_CV20.data)
PM_3pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt_CV20_data.grp)
PM_3pt_CV20.nlme <- nlme(PM_3pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_CV20_covar.nlme <- update(PM_3pt_CV20.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt_CV20.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt_CV20.nlme)[[2]], rep(0,1)))
PM_4pt_CV20.data <- update_data(n_indiv = 10, `CV%` = 20, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt_CV20_data.grp <- groupedData(V~S|ID, PM_4pt_CV20.data)
PM_4pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt_CV20_data.grp)
PM_4pt_CV20.nlme <- nlme(PM_4pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_CV20_covar.nlme <- update(PM_4pt_CV20.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt_CV20.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt_CV20.nlme)[[2]], rep(0,1)))plot(PM_covar.nlme, main = "PM Rich (16pt, 0% CV)")plot(PM_covar.nlme %>% update(weights = varPower()),
main = "PM Rich (16pt, 0% CV, weighted)")plot(PM_3pt_covar.nlme, main = "PM Sparse (3pt, 0% CV)")plot(PM_3pt_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 0% CV, weighted)")plot(PM_4pt_covar.nlme, main = "PM Sparse (4pt, 0% CV)")plot(PM_4pt_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 0% CV, weighted)")plot(PM_CV5_covar.nlme, main = "PM Rich (16pt, 5% CV)")plot(PM_CV5_covar.nlme %>% update(weights = varPower()),
main = "Rich (16pt, 5% CV, weighted)")plot(PM_3pt_CV5_covar.nlme, main = "PM Sparse (3pt, 5% CV)")plot(PM_3pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 5% CV, weighted)")plot(PM_4pt_CV5_covar.nlme, main = "PM Sparse (4pt, 5% CV)")plot(PM_4pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 5% CV, weighted)")plot(PM_CV10_covar.nlme, main = "PM Rich (16pt, 10% CV)")plot(PM_CV10_covar.nlme %>% update(weights = varPower()),
main = "PM Rich (16pt, 10% CV, weighted)")plot(PM_3pt_CV10_covar.nlme, main = "PM Sparse (3pt, 10% CV)")plot(PM_3pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 10% CV, weighted)")plot(PM_4pt_CV10_covar.nlme, main = "PM Sparse (4pt, 10% CV)")plot(PM_4pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 10% CV, weighted)")plot(PM_CV20_covar.nlme, main = "PM Rich (16pt, 20% CV)")plot(PM_CV20_covar.nlme %>% update(weights = varPower()),
main = "PM Rich (16pt, 20% CV, weighted)")plot(PM_3pt_CV20_covar.nlme, main = "PM Sparse (3pt, 20% CV)")plot(PM_3pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 20% CV, weighted)")plot(PM_4pt_CV20_covar.nlme, main = "PM Sparse (4pt, 20% CV)")plot(PM_4pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 20% CV, weighted)")PM_3pt_covar.nlme <- PM_3pt_covar.nlme %>% update(weights = varPower())
PM_4pt_covar.nlme <- PM_4pt_covar.nlme %>% update(weights = varPower())
PM_CV5_covar.nlme <- PM_CV5_covar.nlme %>% update(weights = varPower())
PM_3pt_CV5_covar.nlme <- PM_3pt_CV5_covar.nlme %>% update(weights = varPower())
PM_4pt_CV5_covar.nlme <- PM_4pt_CV5_covar.nlme %>% update(weights = varPower())
PM_CV10_covar.nlme <- PM_CV10_covar.nlme %>% update(weights = varPower())
PM_3pt_CV10_covar.nlme <- PM_3pt_CV10_covar.nlme %>% update(weights = varPower())
PM_4pt_CV10_covar.nlme <- PM_4pt_CV10_covar.nlme %>% update(weights = varPower())
PM_CV20_covar.nlme <- PM_CV20_covar.nlme %>% update(weights = varPower())
PM_3pt_CV20_covar.nlme <- PM_3pt_CV20_covar.nlme %>% update(weights = varPower())
PM_4pt_CV20_covar.nlme <- PM_4pt_CV20_covar.nlme %>% update(weights = varPower())# Summary Data Table ------------
PM_covaref_CV.table <- rbind(
fixef(PM_covar.nlme),
fixef(PM_3pt_covar.nlme),
fixef(PM_4pt_covar.nlme),
fixef(PM_CV5_covar.nlme),
fixef(PM_3pt_CV5_covar.nlme),
fixef(PM_4pt_CV5_covar.nlme),
fixef(PM_CV10_covar.nlme),
fixef(PM_3pt_CV10_covar.nlme),
fixef(PM_4pt_CV10_covar.nlme),
fixef(PM_CV20_covar.nlme),
fixef(PM_3pt_CV20_covar.nlme),
fixef(PM_4pt_CV20_covar.nlme)) %>%
as_tibble() %>%
mutate(`Covariate Model` = c("Rich","Sparse (3pt)","Sparse (4pt)",
"Rich (5% CV)","Sparse (5% CV, 3pt)","Sparse (5% CV, 4pt)",
"Rich (10% CV)","Sparse (10% CV, 3pt)","Sparse (10% CV, 4pt)",
"Rich (20% CV)","Sparse (20% CV, 3pt)","Sparse (20% CV, 4pt)")) %>%
pivot_longer(`Vmax.(Intercept)`:`Km.Diplotype4/5`,
names_to = "Parameter",
values_to = "Value") %>%
pivot_wider(names_from = `Covariate Model`, values_from = Value) %>%
separate(Parameter, sep = "\\.", into = c("Parameter","Diplotype"), remove = T) %>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (4/41)"),
across(where(is.numeric), round, 2))
# Complete Covariate Effects Table --------------------
PM_covar_Vmax.table <- PM_covaref_CV.table %>%filter(Parameter == "Vmax")
PM_covar_Km.table <- PM_covaref_CV.table %>% filter(Parameter == "Km")
# Vmax Estimates -------------------------------
PM_final_vmax.est <- PM_covar_Vmax.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (4/41)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (4/41)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (4/41)" = "Diplotype4/41")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Km), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Vmax)/Vmax)*100, digits = 2),
Vmax = round(Vmax, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Vmax, Predicted, `Error (%)`) %>%
rename("Predicted (Vmax)" = "Predicted") %>%
arrange(Condition)
# Km Estimates -----------------------------------
PM_final_km.est <- PM_covar_Km.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (4/41)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (4/41)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (4/41)" = "Diplotype4/41")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Vmax), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Km)/Km)*100, digits = 2),
Km = round(Km, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Km, Predicted, `Error (%)`) %>%
rename("Predicted (Km)" = "Predicted") %>%
arrange(Condition)all_km.est <- rbind(final_km.est, PM_final_km.est)
all_vmax.est <- rbind(final_vmax.est, PM_final_vmax.est)Vmax.plot <- ggplot(data = all_vmax.est, aes(x = Vmax, y = `Predicted (Vmax)`))+
geom_point(size = 3, alpha = 0.85, aes(fill = Diplotype), shape = 21)+
geom_abline(slope = 1, intercept = 0, color = "red")+
geom_abline(slope = 1, intercept = 0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.2, color = "red", linetype = "dashed")+
scale_y_log10()+
scale_x_log10()+
facet_wrap(.~Condition, ncol = 4, scales = "free")+
xlab("Actual (Vmax)")+
theme_bw(base_size = 12)+
ggeasy::easy_move_legend("top")+
scale_fill_viridis_d()
Vmax.plotKm.plot <- ggplot(data = all_km.est, aes(x = Km, y = `Predicted (Km)`))+
geom_point(size = 3, alpha = 0.85, aes(fill = Diplotype), shape = 21)+
geom_abline(slope = 1, intercept = 0, color = "red")+
geom_abline(slope = 1, intercept = 0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.2, color = "red", linetype = "dashed")+
scale_y_log10()+
scale_x_log10()+
facet_wrap(.~Condition, ncol = 4, scales = "free")+
xlab("Actual (Km)")+
theme_bw(base_size = 12)+
ggeasy::easy_move_legend("top")+
scale_fill_viridis_d()
Km.plot# Function to grab t-Table output from NLME model objects -----------
get_tTable <- function(data.nlme, condition, PM = FALSE){
label <- if_else(PM == TRUE, "4/41", "1/1")
summary(data.nlme)$tTable %>%
as_tibble() %>%
mutate(Parameter = row.names(summary(data.nlme)$tTable)) %>%
relocate(Parameter) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = paste0("Diplotype",label))) %>%
group_by(Parameter) %>%
mutate(Value = if_else(Diplotype == paste0("Diplotype",label),
Value,
Value+Value[Diplotype == paste0("Diplotype",label)])) %>%
separate(Diplotype, remove = T, sep = "Diplotype", into = c(".", "Diplotype")) %>%
select(Parameter, Diplotype, Value, Std.Error) %>%
mutate(`RSE(%)` = round((Std.Error/Value)*100, digits = 2),
Condition = condition)
}# Standard Error Tables
## EM and UM Population -------
rich_CV0_covar.tTable <- get_tTable(rich_covar.nlme, "Rich (0% CV)")
sparse3pt_CV0_covar.tTable <- get_tTable(sparse3pt_covar.nlme, "Sparse (0% CV, 3pt)")
sparse4pt_CV0_covar.tTable <- get_tTable(sparse4pt_covar.nlme, "Sparse (0% CV, 4pt)")
rich_CV5_covar.tTable <- get_tTable(rich_CV5_covar.nlme, "Rich (5% CV)")
sparse3pt_CV5_covar.tTable <- get_tTable(sparse3pt_CV5_covar.nlme, "Sparse (5% CV, 3pt)")
sparse4pt_CV5_covar.tTable <- get_tTable(sparse4pt_CV5_covar.nlme, "Sparse (5% CV, 4pt)")
rich_CV10_covar.tTable <- get_tTable(rich_CV10_covar.nlme, "Rich (10% CV)")
sparse3pt_CV10_covar.tTable <- get_tTable(sparse3pt_CV10_covar.nlme, "Sparse (10% CV, 3pt)")
sparse4pt_CV10_covar.tTable <- get_tTable(sparse4pt_CV10_covar.nlme, "Sparse (10% CV, 4pt)")
rich_CV20_covar.tTable <- get_tTable(rich_CV20_covar.nlme, "Rich (20% CV)")
sparse3pt_CV20_covar.tTable <- get_tTable(sparse3pt_CV20_covar.nlme, "Sparse (20% CV, 3pt)")
sparse4pt_CV20_covar.tTable <- get_tTable(sparse4pt_CV20_covar.nlme, "Sparse (20% CV, 4pt)")
## PM Population -------
PM_CV0_covar.tTable <- get_tTable(PM_covar.nlme, "Rich (0% CV)", T)
PM_3pt_CV0_covar.tTable <- get_tTable(PM_3pt_covar.nlme, "Sparse (0% CV, 3pt)", T)
PM_4pt_CV0_covar.tTable <- get_tTable(PM_4pt_covar.nlme, "Sparse (0% CV, 4pt)", T)
PM_CV5_covar.tTable <- get_tTable(PM_CV5_covar.nlme, "Rich (5% CV)", T)
PM_3pt_CV5_covar.tTable <- get_tTable(PM_3pt_CV5_covar.nlme, "Sparse (5% CV, 3pt)", T)
PM_4pt_CV5_covar.tTable <- get_tTable(PM_4pt_CV5_covar.nlme, "Sparse (5% CV, 4pt)", T)
PM_CV10_covar.tTable <- get_tTable(PM_CV10_covar.nlme, "Rich (10% CV)", T)
PM_3pt_CV10_covar.tTable <- get_tTable(PM_3pt_CV10_covar.nlme, "Sparse (10% CV, 3pt)", T)
PM_4pt_CV10_covar.tTable <- get_tTable(PM_4pt_CV10_covar.nlme, "Sparse (10% CV, 4pt)", T)
PM_CV20_covar.tTable <- get_tTable(PM_CV20_covar.nlme, "Rich (20% CV)", T)
PM_3pt_CV20_covar.tTable <- get_tTable(PM_3pt_CV20_covar.nlme, "Sparse (20% CV, 3pt)", T)
PM_4pt_CV20_covar.tTable <- get_tTable(PM_4pt_CV20_covar.nlme, "Sparse (20% CV, 4pt)", T)## Combined tTable Outpout -----------------
complete_CV_covar.tTable <- rbind(
# EM and UM Population Estimates -------
rich_CV5_covar.tTable,
sparse3pt_CV5_covar.tTable,
sparse4pt_CV5_covar.tTable,
rich_CV10_covar.tTable,
sparse3pt_CV10_covar.tTable,
sparse4pt_CV10_covar.tTable,
rich_CV20_covar.tTable,
sparse3pt_CV20_covar.tTable,
sparse4pt_CV20_covar.tTable,
#PM and IM Population Estimates -----
PM_CV5_covar.tTable,
PM_3pt_CV5_covar.tTable,
PM_4pt_CV5_covar.tTable,
PM_CV10_covar.tTable,
PM_3pt_CV10_covar.tTable,
PM_4pt_CV10_covar.tTable,
PM_CV20_covar.tTable,
PM_3pt_CV20_covar.tTable,
PM_4pt_CV20_covar.tTable)# Confidence Interval Table ---------
## EM and UM Population Intervals
rich_CV5_covar.intervals <- intervals(rich_CV5_covar.nlme, which = "fixed")
sparse3pt_CV5_covar.intervals <- intervals(sparse3pt_CV5_covar.nlme, which = "fixed")
sparse4pt_CV5_covar.intervals <- intervals(sparse4pt_CV5_covar.nlme, which = "fixed")
rich_CV10_covar.intervals <- intervals(rich_CV10_covar.nlme, which = "fixed")
sparse3pt_CV10_covar.intervals <- intervals(sparse3pt_CV10_covar.nlme, which = "fixed")
sparse4pt_CV10_covar.intervals <- intervals(sparse4pt_CV10_covar.nlme, which = "fixed")
rich_CV20_covar.intervals <- intervals(rich_CV20_covar.nlme, which = "fixed")
sparse3pt_CV20_covar.intervals <- intervals(sparse3pt_CV20_covar.nlme, which = "fixed")
sparse4pt_CV20_covar.intervals <- intervals(sparse4pt_CV20_covar.nlme, which = "fixed")
## PM and IM Population Intervals
PM_covar.intervals <- intervals(PM_covar.nlme, which = "fixed")
PM_3pt_covar.intervals <- intervals(PM_3pt_covar.nlme, which = "fixed")
PM_4pt_covar.intervals <- intervals(PM_4pt_covar.nlme, which = "fixed")
PM_CV5_covar.intervals <- intervals(PM_CV5_covar.nlme, which = "fixed")
PM_3pt_CV5_covar.intervals <- intervals(PM_3pt_CV5_covar.nlme, which = "fixed")
PM_4pt_CV5_covar.intervals <- intervals(PM_4pt_CV5_covar.nlme, which = "fixed")
PM_CV10_covar.intervals <- intervals(PM_CV10_covar.nlme, which = "fixed")
PM_3pt_CV10_covar.intervals <- intervals(PM_3pt_CV10_covar.nlme, which = "fixed")
PM_4pt_CV10_covar.intervals <- intervals(PM_4pt_CV10_covar.nlme, which = "fixed")
PM_CV20_covar.intervals <- intervals(PM_CV20_covar.nlme, which = "fixed")
PM_3pt_CV20_covar.intervals <- intervals(PM_3pt_CV20_covar.nlme, which = "fixed")
PM_4pt_CV20_covar.intervals <- intervals(PM_4pt_CV20_covar.nlme, which = "fixed")# Function to extract and tidy confidence interval estimates ---------
extract_CI3 <- function(int.data, name, PM = FALSE){
label <- if_else(PM == TRUE, "4/41", "1/1")
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(Parameter = rownames(int.data$fixed)) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = paste0("Diplotype",label))) %>%
group_by(Parameter) %>%
mutate(lower = if_else(Diplotype != paste0("Diplotype",label),
lower + est.[Diplotype == paste0("Diplotype",label)],
lower),
upper = if_else(Diplotype != paste0("Diplotype",label),
upper + est.[Diplotype == paste0("Diplotype",label)],
upper),
est. = if_else(Diplotype != paste0("Diplotype",label),
est. + est.[Diplotype == paste0("Diplotype",label)],
est.),
across(where(is.numeric), round, 2),
`CI (95%)` = paste0("(",lower,", ",upper,")")) %>%
ungroup() %>%
separate(Diplotype, remove = T, sep = "Diplotype", into = c(".", "Diplotype"))%>%
select(Parameter, Diplotype, est., `CI (95%)`)
colnames(temp) <- c("Parameter", "Diplotype",
paste(name, " Estimate"),
paste(name, " CI (95%)"))
temp
}# Combined Tables
CV0_covar_ci.table <- cbind(
extract_CI3(rich_covar.intervals, name = "Rich (0% CV)"),
rich_CV0_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_covar.intervals, name = "Sparse 3pt (0% CV)")[,3:4],
sparse3pt_CV0_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_covar.intervals, name = "Sparse 4pt (0% CV)")[,3:4],
sparse4pt_CV0_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
CV5_covar_ci.table <- cbind(
extract_CI3(rich_CV5_covar.intervals, name = "Rich (5% CV)"),
rich_CV5_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_CV5_covar.intervals, name = "Sparse 3pt (5% CV)")[,3:4],
sparse3pt_CV5_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_CV5_covar.intervals, name = "Sparse 4pt (5% CV)")[,3:4],
sparse4pt_CV5_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
CV10_covar_ci.table <- cbind(
extract_CI3(rich_CV10_covar.intervals, name = "Rich (10% CV)"),
rich_CV10_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_CV10_covar.intervals, name = "Sparse 3pt (10% CV)")[,3:4],
sparse3pt_CV10_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_CV10_covar.intervals, name = "Sparse 4pt (10% CV)")[,3:4],
sparse4pt_CV10_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
CV20_covar_ci.table <- cbind(
extract_CI3(rich_CV20_covar.intervals, name = "Rich (20% CV)"),
rich_CV20_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_CV20_covar.intervals, name = "Sparse 3pt (20% CV)")[,3:4],
sparse3pt_CV20_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_CV20_covar.intervals, name = "Sparse 4pt (20% CV)")[,3:4],
sparse4pt_CV20_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV0_covar_ci.table <- cbind(
extract_CI3(PM_covar.intervals, name = "Rich (0% CV)", T),
PM_CV0_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_covar.intervals, name = "Sparse 3pt (0% CV)", T)[,3:4],
PM_3pt_CV0_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_covar.intervals, name = "Sparse 4pt (0% CV)", T)[,3:4],
PM_4pt_CV0_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV5_covar_ci.table <- cbind(
extract_CI3(PM_CV5_covar.intervals, name = "Rich (5% CV)", T),
PM_CV5_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_CV5_covar.intervals, name = "Sparse 3pt (5% CV)", T)[,3:4],
PM_3pt_CV5_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_CV5_covar.intervals, name = "Sparse 4pt (5% CV)", T)[,3:4],
PM_4pt_CV5_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV10_covar_ci.table <- cbind(
extract_CI3(PM_CV10_covar.intervals, name = "Rich (10% CV)", T),
PM_CV10_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_CV10_covar.intervals, name = "Sparse 3pt (10% CV)", T)[,3:4],
PM_3pt_CV10_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_CV10_covar.intervals, name = "Sparse 4pt (10% CV)", T)[,3:4],
PM_4pt_CV10_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV20_covar_ci.table <- cbind(
extract_CI3(PM_CV20_covar.intervals, name = "Rich (20% CV)", T),
PM_CV20_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_CV20_covar.intervals, name = "Sparse 3pt (20% CV)", T)[,3:4],
PM_3pt_CV20_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_CV20_covar.intervals, name = "Sparse 4pt (20% CV)", T)[,3:4],
PM_4pt_CV20_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)colnames(CV0_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(CV5_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(CV10_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(CV20_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV0_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV5_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV10_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV20_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")CV0_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 0% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich" = 3,
"Sparse (3pt)" = 3,
"Sparse (4pt)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich
|
Sparse (3pt)
|
Sparse (4pt)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 367.01 | (319.69, 414.33) | 6.64 | 367.01 | (319.34, 414.67) | 6.79 | 367.01 | (319.51, 414.51) | 6.73 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 486.93 | (420.01, 553.85) | 7.08 | 486.93 | (419.52, 554.34) | 7.24 | 486.93 | (419.75, 554.11) | 7.18 |
| Vmax | 1/2x2 | 1410.00 | 1459.75 | (1392.84, 1526.67) | 2.36 | 1459.75 | (1392.34, 1527.16) | 2.42 | 1459.75 | (1392.57, 1526.93) | 2.39 |
| Vmax | 1/4 | 274.00 | 283.67 | (216.75, 350.58) | 12.14 | 283.67 | (216.26, 351.08) | 12.43 | 283.67 | (216.49, 350.85) | 12.32 |
| Vmax | 2/2 | 252.00 | 260.89 | (193.98, 327.81) | 13.21 | 260.89 | (193.48, 328.3) | 13.52 | 260.89 | (193.71, 328.07) | 13.40 |
| Vmax | 2/3 | 251.00 | 259.86 | (192.94, 326.77) | 13.26 | 259.86 | (192.45, 327.27) | 13.57 | 259.86 | (192.68, 327.04) | 13.45 |
| Vmax | 2/4 | 95.90 | 99.28 | (32.37, 166.2) | 34.70 | 99.28 | (31.87, 166.69) | 35.52 | 99.28 | (32.1, 166.46) | 35.20 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.50 | (1.33, 1.67) | 5.92 | 1.50 | (1.33, 1.67) | 6.06 | 1.50 | (1.33, 1.67) | 6.00 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.44 | (1.2, 1.68) | 8.73 | 1.44 | (1.19, 1.68) | 8.93 | 1.44 | (1.19, 1.68) | 8.85 |
| Km | 1/2x2 | 1.82 | 1.89 | (1.65, 2.13) | 6.65 | 1.89 | (1.64, 2.14) | 6.80 | 1.89 | (1.64, 2.13) | 6.74 |
| Km | 1/4 | 0.92 | 0.95 | (0.71, 1.2) | 13.18 | 0.95 | (0.71, 1.2) | 13.49 | 0.95 | (0.71, 1.2) | 13.37 |
| Km | 2/2 | 4.59 | 4.75 | (4.51, 5) | 2.64 | 4.75 | (4.51, 5) | 2.71 | 4.75 | (4.51, 5) | 2.68 |
| Km | 2/3 | 1.32 | 1.37 | (1.12, 1.61) | 9.19 | 1.37 | (1.12, 1.61) | 9.41 | 1.37 | (1.12, 1.61) | 9.32 |
| Km | 2/4 | 1.68 | 1.74 | (1.49, 1.98) | 7.24 | 1.74 | (1.49, 1.98) | 7.41 | 1.74 | (1.49, 1.98) | 7.34 |
CV5_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 5% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (5% CV)" = 3,
"Sparse (3pt, 5% CV)" = 3,
"Sparse (4pt, 5% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich (5% CV)
|
Sparse (3pt, 5% CV)
|
Sparse (4pt, 5% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 362.01 | (326.7, 397.31) | 5.02 | 364.79 | (346.49, 383.1) | 2.63 | 363.04 | (343.06, 383.02) | 2.86 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 487.84 | (437.63, 538.06) | 5.30 | 482.00 | (452.75, 511.25) | 3.17 | 481.46 | (451.75, 511.16) | 3.21 |
| Vmax | 1/2x2 | 1410.00 | 1461.41 | (1405.98, 1516.83) | 1.95 | 1463.40 | (1387.76, 1539.04) | 2.70 | 1460.85 | (1405.71, 1515.98) | 1.96 |
| Vmax | 1/4 | 274.00 | 283.50 | (233.72, 333.28) | 9.04 | 278.28 | (254.56, 302.01) | 4.46 | 278.34 | (250.88, 305.8) | 5.13 |
| Vmax | 2/2 | 252.00 | 261.57 | (211.72, 311.41) | 9.81 | 264.06 | (239.64, 288.48) | 4.84 | 264.62 | (236.97, 292.26) | 5.44 |
| Vmax | 2/3 | 251.00 | 261.41 | (211.64, 311.17) | 9.80 | 258.97 | (235.43, 282.51) | 4.75 | 259.09 | (231.73, 286.45) | 5.49 |
| Vmax | 2/4 | 95.90 | 97.83 | (48.22, 147.45) | 26.11 | 98.79 | (77.14, 120.44) | 11.46 | 98.20 | (71.53, 124.86) | 14.13 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.47 | (1.34, 1.59) | 4.26 | 1.49 | (1.38, 1.59) | 3.74 | 1.47 | (1.36, 1.58) | 3.87 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.46 | (1.29, 1.63) | 6.03 | 1.46 | (1.31, 1.62) | 5.41 | 1.47 | (1.31, 1.62) | 5.54 |
| Km | 1/2x2 | 1.82 | 1.88 | (1.71, 2.06) | 4.81 | 1.89 | (1.71, 2.08) | 5.18 | 1.90 | (1.72, 2.07) | 4.87 |
| Km | 1/4 | 0.92 | 0.97 | (0.8, 1.13) | 8.93 | 0.96 | (0.83, 1.09) | 6.97 | 0.94 | (0.79, 1.08) | 7.96 |
| Km | 2/2 | 4.59 | 4.83 | (4.61, 5.05) | 2.37 | 4.99 | (4.65, 5.32) | 3.51 | 4.99 | (4.73, 5.25) | 2.71 |
| Km | 2/3 | 1.32 | 1.39 | (1.22, 1.56) | 6.33 | 1.37 | (1.22, 1.51) | 5.48 | 1.37 | (1.22, 1.52) | 5.72 |
| Km | 2/4 | 1.68 | 1.68 | (1.5, 1.85) | 5.31 | 1.67 | (1.52, 1.81) | 4.61 | 1.65 | (1.5, 1.8) | 4.73 |
CV10_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 10% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (10% CV)" = 3,
"Sparse (3pt, 10% CV)" = 3,
"Sparse (4pt, 10% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich (10% CV)
|
Sparse (3pt, 10% CV)
|
Sparse (4pt, 10% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 355.96 | (336.28, 375.64) | 2.85 | 365.54 | (341.6, 389.48) | 3.43 | 361.21 | (339.36, 383.05) | 3.15 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 487.46 | (457.38, 517.54) | 3.18 | 481.34 | (442.9, 519.79) | 4.18 | 479.61 | (445.64, 513.58) | 3.68 |
| Vmax | 1/2x2 | 1410.00 | 1479.15 | (1420.74, 1537.56) | 2.03 | 1476.20 | (1380.61, 1571.8) | 3.39 | 1474.03 | (1400.08, 1547.98) | 2.61 |
| Vmax | 1/4 | 274.00 | 282.04 | (255.33, 308.75) | 4.88 | 273.99 | (243.4, 304.57) | 5.84 | 275.03 | (246.08, 303.98) | 5.48 |
| Vmax | 2/2 | 252.00 | 263.89 | (236.56, 291.22) | 5.33 | 267.54 | (235.45, 299.63) | 6.27 | 268.41 | (238.59, 298.23) | 5.78 |
| Vmax | 2/3 | 251.00 | 262.09 | (235.47, 288.71) | 5.23 | 259.86 | (229.39, 290.33) | 6.13 | 259.08 | (230.28, 287.89) | 5.78 |
| Vmax | 2/4 | 95.90 | 96.22 | (70.87, 121.57) | 13.57 | 98.98 | (71.55, 126.42) | 14.50 | 97.58 | (70.7, 124.47) | 14.33 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.43 | (1.33, 1.53) | 3.63 | 1.48 | (1.32, 1.65) | 5.85 | 1.45 | (1.3, 1.61) | 5.43 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.49 | (1.34, 1.63) | 5.03 | 1.51 | (1.28, 1.75) | 8.22 | 1.51 | (1.29, 1.72) | 7.54 |
| Km | 1/2x2 | 1.82 | 1.91 | (1.74, 2.07) | 4.51 | 1.92 | (1.64, 2.2) | 7.51 | 1.92 | (1.67, 2.17) | 6.74 |
| Km | 1/4 | 0.92 | 0.97 | (0.85, 1.1) | 6.55 | 0.96 | (0.76, 1.17) | 10.98 | 0.96 | (0.78, 1.15) | 10.12 |
| Km | 2/2 | 4.59 | 4.97 | (4.63, 5.31) | 3.54 | 5.19 | (4.64, 5.74) | 5.54 | 5.19 | (4.7, 5.69) | 4.98 |
| Km | 2/3 | 1.32 | 1.41 | (1.26, 1.55) | 5.17 | 1.38 | (1.15, 1.6) | 8.61 | 1.37 | (1.16, 1.57) | 7.95 |
| Km | 2/4 | 1.68 | 1.61 | (1.46, 1.76) | 4.76 | 1.62 | (1.38, 1.86) | 7.70 | 1.59 | (1.37, 1.8) | 7.19 |
CV20_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 20% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (20% CV)" = 3,
"Sparse (3pt, 20% CV)" = 3,
"Sparse (4pt, 20% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich (20% CV)
|
Sparse (3pt, 20% CV)
|
Sparse (4pt, 20% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 345.88 | (319.99, 371.77) | 3.85 | 366.90 | (329.06, 404.73) | 5.39 | 358.47 | (326.75, 390.19) | 4.60 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 491.31 | (448.44, 534.19) | 4.49 | 484.18 | (421.45, 546.91) | 6.78 | 477.51 | (426, 529.02) | 5.61 |
| Vmax | 1/2x2 | 1410.00 | 1510.05 | (1408.03, 1612.07) | 3.48 | 1501.07 | (1339.29, 1662.86) | 5.64 | 1497.51 | (1373.43, 1621.59) | 4.31 |
| Vmax | 1/4 | 274.00 | 280.85 | (246.74, 314.97) | 6.25 | 265.12 | (218.62, 311.62) | 9.18 | 267.56 | (227.18, 307.95) | 7.85 |
| Vmax | 2/2 | 252.00 | 266.92 | (230.87, 302.98) | 6.95 | 275.51 | (224.49, 326.53) | 9.69 | 275.49 | (232.06, 318.91) | 8.20 |
| Vmax | 2/3 | 251.00 | 264.38 | (230.39, 298.38) | 6.62 | 262.08 | (215.29, 308.86) | 9.34 | 259.33 | (218.97, 299.7) | 8.10 |
| Vmax | 2/4 | 95.90 | 93.65 | (63.33, 123.97) | 16.67 | 98.20 | (58.18, 138.22) | 21.32 | 95.77 | (60.19, 131.34) | 19.32 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.37 | (1.18, 1.55) | 6.92 | 1.48 | (1.15, 1.8) | 11.68 | 1.43 | (1.14, 1.72) | 10.53 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.56 | (1.28, 1.83) | 9.13 | 1.64 | (1.16, 2.13) | 15.47 | 1.61 | (1.18, 2.04) | 13.84 |
| Km | 1/2x2 | 1.82 | 1.94 | (1.63, 2.25) | 8.25 | 1.98 | (1.45, 2.51) | 14.06 | 1.98 | (1.51, 2.45) | 12.32 |
| Km | 1/4 | 0.92 | 1.00 | (0.77, 1.23) | 11.91 | 0.96 | (0.54, 1.37) | 22.61 | 0.98 | (0.62, 1.34) | 19.16 |
| Km | 2/2 | 4.59 | 5.20 | (4.52, 5.88) | 6.71 | 5.62 | (4.48, 6.76) | 10.62 | 5.58 | (4.53, 6.64) | 9.82 |
| Km | 2/3 | 1.32 | 1.45 | (1.18, 1.72) | 9.47 | 1.41 | (0.95, 1.86) | 17.00 | 1.38 | (0.98, 1.78) | 15.24 |
| Km | 2/4 | 1.68 | 1.51 | (1.24, 1.79) | 9.29 | 1.53 | (1.06, 2) | 16.19 | 1.46 | (1.04, 1.88) | 14.96 |
PM_CV0_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 0% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich" = 3,
"Sparse (3pt)" = 3,
"Sparse (4pt)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich
|
Sparse (3pt)
|
Sparse (4pt)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 31.99 | (30.14, 33.84) | 2.96 | 31.99 | (30.28, 33.69) | 2.72 | 31.94 | (30.2, 33.67) | 2.79 |
| Vmax | 4/5 | 12.32 | 12.75 | (10.14, 15.37) | 10.48 | 12.74 | (10.34, 15.15) | 9.65 | 12.43 | (10.64, 14.23) | 7.40 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.54 | (1.69, 9.39) | 35.57 | 5.61 | (5.3, 5.92) | 2.83 | 5.54 | (5.07, 6.02) | 4.38 |
| Km | 4/5 | 69.15 | 71.60 | (66.15, 77.05) | 3.89 | 71.96 | (67.76, 76.15) | 2.98 | 68.35 | (65.66, 71.04) | 2.02 |
PM_CV5_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 5% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (5% CV)" = 3,
"Sparse (3pt, 5% CV)" = 3,
"Sparse (4pt, 5% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich (5% CV)
|
Sparse (3pt, 5% CV)
|
Sparse (4pt, 5% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 31.70 | (30.04, 33.35) | 2.67 | 31.88 | (30.1, 33.65) | 2.85 | 31.62 | (29.89, 33.36) | 2.81 |
| Vmax | 4/5 | 12.32 | 12.60 | (10.27, 14.92) | 9.44 | 12.47 | (10.11, 14.84) | 9.70 | 12.56 | (10.72, 14.4) | 7.51 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.56 | (5.17, 5.96) | 3.65 | 6.19 | (5.07, 7.31) | 9.23 | 5.57 | (5.03, 6.11) | 4.97 |
| Km | 4/5 | 69.15 | 69.17 | (67.09, 71.25) | 1.54 | 67.23 | (61.18, 73.28) | 4.60 | 69.26 | (65.67, 72.85) | 2.66 |
PM_CV10_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 10% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (10% CV)" = 3,
"Sparse (3pt, 10% CV)" = 3,
"Sparse (4pt, 10% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich (10% CV)
|
Sparse (3pt, 10% CV)
|
Sparse (4pt, 10% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 31.37 | (29.86, 32.87) | 2.45 | 32.02 | (29.53, 34.52) | 3.98 | 31.40 | (29.23, 33.57) | 3.55 |
| Vmax | 4/5 | 12.32 | 12.50 | (10.44, 14.56) | 8.44 | 12.25 | (9.58, 14.92) | 11.12 | 12.46 | (10.1, 14.81) | 9.68 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.60 | (5.2, 5.99) | 3.60 | 6.82 | (4.18, 9.46) | 19.74 | 5.66 | (4.87, 6.45) | 7.16 |
| Km | 4/5 | 69.15 | 67.53 | (64.29, 70.77) | 2.46 | 63.57 | (54.88, 72.26) | 6.99 | 68.20 | (62.5, 73.89) | 4.28 |
PM_CV20_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 20% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (20% CV)" = 3,
"Sparse (3pt, 20% CV)" = 3,
"Sparse (4pt, 20% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich (20% CV)
|
Sparse (3pt, 20% CV)
|
Sparse (4pt, 20% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 30.77 | (29.18, 32.36) | 2.64 | 32.17 | (28.23, 36.11) | 6.26 | 31.02 | (27.58, 34.46) | 5.68 |
| Vmax | 4/5 | 12.32 | 12.33 | (10.37, 14.3) | 8.16 | 11.89 | (7.68, 16.1) | 18.10 | 12.28 | (8.52, 16.03) | 15.67 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.70 | (4.93, 6.48) | 6.96 | 8.33 | (3.77, 12.89) | 27.95 | 5.85 | (4.43, 7.28) | 12.49 |
| Km | 4/5 | 69.15 | 64.82 | (59.1, 70.55) | 4.51 | 57.03 | (42.81, 71.26) | 12.74 | 65.65 | (55.21, 76.09) | 8.15 |
The predictor function sim_nlme() was scripted to
generate tidy data frames of NLME predictions given the specified
substrate concentration range. The inputs for this function are as
followed:
tTable - extracted NLME model tTable estimates (see Final Model Parameter Estimates and Standard Error Datasets section).
label - Specifies the condition being simulated (ex. “Rich 0% CV”).
Group - Specifies the CYP2D6 metabolizer category (ex. “IM & PM”).
Start - Substrate concentration to start simulation (default = 0 uM).
End - Last substrate concentration (default = 100 uM).
by - Simulated concentration increment (default = 0.1)
sim_nlme <- function(tTable, label, Group, Start = 0, End = 100, by = 0.1){
diplotypes_contained <- unique(tTable$Diplotype)
# Simulate Diplotype Rates -----
tTable %>%
select(Parameter, Diplotype, Value) %>%
pivot_wider(names_from = Parameter,
values_from = Value) %>%
expand_grid(S = seq(Start,End,by)) %>%
mutate(V = (Vmax*S)/(Km+S),
Condition = label,
Group = Group)
}update_data_full() is an extended version of the
previous update_data() function; it differs by retaining Km
and Vmax information in the output. It was created as an intermediate
function for the revised update_data2() function (see
below).
update_data_full <-function(`CV%`, type, n_indiv = 10, `CV%.D` = 25, by = 0.1,
start = 0, end = 100, Seed = 23457, Pop.freq = FALSE){
# Sampling Information --------------
# Full Range Set
full_set <- c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)
PM_range <- c(1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, 2000)
# Strategic Sampling Sets ----
sparse3pt_set1 <- c(0.1, 2.5, 50)
sparse3pt_set2 <- c(1, 10, 100)
sparse4pt_set1 <- c(0.1, 2.5, 25, 50)
sparse4pt_set2 <- c(1, 10, 50, 100)
PM_3pt_set1 <- c(10,100,1000)
PM_3pt_set2 <- c(25,250,2000)
PM_4pt_set1 <- c(1,10,100,1000)
PM_4pt_set2 <- c(5,25,250,2000)
PM_range_set1 <- c(1,10,25,50,100,400,1000,2000)
PM_range_set2 <- c(5,15,30,60,90,200,800,1600)
# Rich Sampling Simulation ----------------
Rich <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
# Sparse Sampling Simulation (3 point) ------------
Sparse3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse3pt_set1, S %in% sparse3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
# Sparse Sampling Simulation (4 point) --------------
Sparse4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse4pt_set1, S %in% sparse4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
PM <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
PM_3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_3pt_set1, S %in% PM_3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
PM_4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_4pt_set1, S %in% PM_4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
# Output Selection ----------
switch(type,
"rich" = Rich,
"sparse3pt" = Sparse3pt,
"sparse4pt" = Sparse4pt,
"PM" = PM,
"PM_3pt" = PM_3pt,
"PM_4pt" = PM_4pt)
}Revised vs of update_data() designed to allow user to
include more information (i.e., Condition, and Group) to the simulated
data output.
update_data2 <- function(`CV%`, type, Condition, Group,
n_indiv = 10,
`CV%.D` = 25,
by = 0.1,
start = 0,
end = 100,
Seed = 23457,
Pop.freq = FALSE){
update_data_full(`CV%`, type, n_indiv = n_indiv,
`CV%.D` = `CV%.D`, by = by,
start = start, end = end,
Seed = Seed, Pop.freq = Pop.freq) %>%
mutate(Condition = Condition,
Group = Group) %>%
filter(if_else(Group == "EM & UM",
!Diplotype %in% c("4/41", "4/5"),
Diplotype %in% c("4/41", "4/5")))
}# Population Level Simulated Data Based on Experimental Conditions ---
sensitivity.data <- rbind(sim_nlme(rich_CV0_covar.tTable,"Rich (0% CV)", "EM & UM"),
sim_nlme(rich_CV5_covar.tTable,"Rich (5% CV)", "EM & UM"),
sim_nlme(rich_CV10_covar.tTable,"Rich (10% CV)", "EM & UM"),
sim_nlme(rich_CV20_covar.tTable,"Rich (20% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV0_covar.tTable,"Sparse (3pt, 0% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV5_covar.tTable,"Sparse (3pt, 5% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV10_covar.tTable,"Sparse (3pt, 10% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV20_covar.tTable,"Sparse (3pt, 20% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV0_covar.tTable,"Sparse (4pt, 0% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV5_covar.tTable,"Sparse (4pt, 5% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV10_covar.tTable,"Sparse (4pt, 10% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV20_covar.tTable,"Sparse (4pt, 20% CV)", "EM & UM"),
sim_nlme(PM_CV0_covar.tTable, "Rich (0% CV)", "IM & PM", End = 2000),
sim_nlme(PM_CV5_covar.tTable, "Rich (5% CV)", "IM & PM", End = 2000),
sim_nlme(PM_CV10_covar.tTable, "Rich (10% CV)", "IM & PM", End = 2000),
sim_nlme(PM_CV20_covar.tTable, "Rich (20% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV0_covar.tTable, "Sparse (3pt, 0% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV5_covar.tTable, "Sparse (3pt, 5% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV10_covar.tTable, "Sparse (3pt, 10% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV20_covar.tTable, "Sparse (3pt, 20% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV0_covar.tTable, "Sparse (4pt, 0% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV5_covar.tTable, "Sparse (4pt, 5% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV10_covar.tTable, "Sparse (4pt, 10% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV20_covar.tTable, "Sparse (4pt, 20% CV)", "IM & PM", End = 2000))
sensitivity.data$Condition <- factor(sensitivity.data$Condition,
levels = unique(sensitivity.data$Condition))# Individual Level Simulated Data Based on Experimental Conditions --
sensitivity.pop <- rbind(update_data2(0,"rich", "Rich (0% CV)", "EM & UM"),
update_data2(5,"rich", "Rich (5% CV)", "EM & UM"),
update_data2(10,"rich", "Rich (10% CV)", "EM & UM"),
update_data2(20,"rich", "Rich (20% CV)", "EM & UM"),
update_data2(0,"sparse3pt", "Sparse (3pt, 0% CV)", "EM & UM"),
update_data2(5,"sparse3pt", "Sparse (3pt, 5% CV)", "EM & UM"),
update_data2(10,"sparse3pt", "Sparse (3pt, 10% CV)", "EM & UM"),
update_data2(20,"sparse3pt", "Sparse (3pt, 20% CV)", "EM & UM"),
update_data2(0,"sparse4pt", "Sparse (4pt, 0% CV)", "EM & UM"),
update_data2(5,"sparse4pt", "Sparse (4pt, 5% CV)", "EM & UM"),
update_data2(10,"sparse4pt", "Sparse (4pt, 10% CV)", "EM & UM"),
update_data2(20,"sparse4pt", "Sparse (4pt, 20% CV)", "EM & UM"),
update_data2(0,"PM", "Rich (0% CV)", "IM & PM"),
update_data2(5,"PM", "Rich (5% CV)", "IM & PM"),
update_data2(10,"PM", "Rich (10% CV)", "IM & PM"),
update_data2(20,"PM", "Rich (20% CV)", "IM & PM"),
update_data2(0,"PM_3pt", "Sparse (3pt, 0% CV)", "IM & PM"),
update_data2(5,"PM_3pt", "Sparse (3pt, 5% CV)", "IM & PM"),
update_data2(10,"PM_3pt", "Sparse (3pt, 10% CV)", "IM & PM"),
update_data2(20,"PM_3pt", "Sparse (3pt, 20% CV)", "IM & PM"),
update_data2(0,"PM_4pt", "Sparse (4pt, 0% CV)", "IM & PM"),
update_data2(5,"PM_4pt", "Sparse (4pt, 5% CV)", "IM & PM"),
update_data2(10,"PM_4pt", "Sparse (4pt, 10% CV)", "IM & PM"),
update_data2(20,"PM_4pt", "Sparse (4pt, 20% CV)", "IM & PM"))
sensitivity.pop$Condition <- factor(sensitivity.pop$Condition,
levels = unique(sensitivity.pop$Condition))sensitivity_mid_end.plot <- ggplot(data = sensitivity.data %>%
filter(S <=100,
!Diplotype %in% c("1/2x2", "4/41", "4/5")),
aes(x = S, y = V, group = Diplotype))+
geom_line(data = sensitivity.pop %>%
filter(S<=100, !Diplotype %in% c("1/2x2", "4/41", "4/5")),
color = "grey",
size = 0.5,
aes(x = S, y = V, group = ID))+
geom_point(data = sensitivity.pop %>%
filter(S<=100, !Diplotype %in% c("1/2x2", "4/41", "4/5")),
aes(x = S, y = V, group = ID, color = Diplotype),
alpha = 0.4,
size = 2)+
geom_line(aes(color = Diplotype), size = 1)+
facet_wrap(~Condition, ncol = 4, nrow = 3, scales = "free")+
theme_bw()+
scale_color_viridis_d()+
xlab("Atomoxetine (uM)")+
ylab("4-OH-Atomoxetine Formation\n(pmol/min/mg protein)")+
ggeasy::easy_move_legend(to = "top")
sensitivity_mid_end.plot+
ggeasy::easy_add_legend_title("Genotype")In-silico NLME predicted Michaelis-Menten fits of 4-hydroxy-atomoxetine formation using human liver microsomes by CYP2D6 extensive metabolizers (1/1, 1/2, 1/4, 2/2, 2/3, and 2/4) across all experimental conditions (sparse 3-4pt, 0-20% CV).
# Rapid Metabolizers -------
sensitivity_high_end.plot <- ggplot(data = sensitivity.data %>%
filter(Diplotype == "1/2x2"),
aes(x = S, y = V, group = Diplotype))+
geom_line(data = sensitivity.pop %>%
filter(Diplotype == "1/2x2"),
color = "grey",
aes(x = S, y = V, group = ID),
size = 1)+
geom_point(data = sensitivity.pop %>%
filter(Diplotype == "1/2x2"),
aes(x = S, y = V, group = ID,color = Diplotype),
alpha = 0.4,
size = 2.5)+
geom_line(aes(color = Diplotype), size = 1.65)+
facet_wrap(Group~Condition, ncol = 4, nrow = 3, scales = "free")+
sjPlot::theme_sjplot()+
scale_color_viridis_d(option = "B", begin = 0)+
xlab("Atomoxetine (uM)")+
ylab("4-OH-Atomoxetine Formation\n(pmol/min/mg protein)")+
ggeasy::easy_move_legend(to = "top")
sensitivity_high_end.plot +
ggeasy::easy_add_legend_title("Genotype")In-silico NLME predicted Michaelis-Menten fits of 4-hydroxy-atomoxetine formation using human liver microsomes by CYP2D6 ultra-rapid metabolizers (1/2X2) across all experimental conditions (sparse 3-4pt, 0-20% CV).
# Poor Metabolizers ----
sensitivity_low_end.plot <- ggplot(data = sensitivity.data %>%
filter(Diplotype %in% c("4/41", "4/5")),
aes(x = S, y = V, group = Diplotype))+
geom_line(data = sensitivity.pop %>%
filter(Diplotype %in% c("4/41", "4/5")),
color = "grey",
size = 1,
aes(x = S, y = V, group = ID))+
geom_line(aes(color = Diplotype), size = 1.65)+
geom_point(data = sensitivity.pop %>%
filter(Diplotype %in% c("4/41", "4/5")),
aes(x = S, y = V, group = ID,color = Diplotype),
alpha = 0.4,
size = 2.5)+
facet_wrap(Group~Condition, ncol = 4, nrow = 3, scales = "free")+
sjPlot::theme_sjplot()+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.5)+
# scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
# labels = trans_format("log10", math_format(10^.x)))+
xlab("Atomoxetine (uM)")+
ylab("4-OH-Atomoxetine Formation\n(pmol/min/mg protein)")+
ggeasy::easy_move_legend(to = "top")
sensitivity_low_end.plot+
ggeasy::easy_add_legend_title("Genotype")In-silico NLME predicted Michaelis-Menten fits of 4-hydroxy-atomoxetine formation using human liver microsomes by CYP2D6 intermediat and poor metabolizers (4/41, and 4/5) across all experimental conditions (sparse 3-4pt, 0-20% CV).
# Data frame with intrinsic clearance estimates for all conditions ----
sensitivity.CL <- sensitivity.data %>%
filter(S == 0) %>%
mutate(CLint = round(Vmax/Km,digits = 2)) %>%
select(-S, -V) %>%
left_join(diplotype.data %>%
mutate(CLint_actual = round(Vmax/Km, digit = 2)) %>%
select(Diplotype, CLint_actual),by = "Diplotype") %>%
mutate(CLint_ratio = CLint/CLint_actual)%>%
group_by(Diplotype) %>%
mutate(standard_resid = (CLint-CLint_actual)/sqrt(CLint_actual),
Obs = seq_along(standard_resid),
`Simulated Error` = if_else(Condition %in% c("Rich (0% CV)",
"Sparse (3pt, 0% CV)",
"Sparse (4pt, 0% CV)"),
"Control\n(0% CV)",
"Variable Error\n(5-20% CV)")) %>%
ungroup()
CL.plot <- ggplot(data = sensitivity.CL, aes(x = CLint_actual, y = CLint))+
geom_point(size = 3.5, alpha = 0.65, aes(fill = Condition), shape = 21)+
geom_abline(slope = 1, intercept = 0, color = "red")+
geom_abline(slope = 1, intercept = 0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = 0.15, color = "blue", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.15, color = "blue", linetype = "dashed")+
scale_y_log10()+
scale_x_log10()+
# facet_wrap(.~Condition, ncol = 4, scales = "free")+
xlab("Intrinsic Clearance (Actual)")+
ylab("Intrinsic Clearance (Predicted)")+
theme_bw(base_size = 14)+
ggeasy::easy_move_legend("right")+
scale_fill_viridis_d()+
# coord_cartesian(clip = "off")+
facet_wrap(~Group, ncol = 1, scales = "free")CL.plot CL_plot2.1 <- ggplot(data = sensitivity.CL,
aes(y = standard_resid, x = Diplotype))+
geom_hline(yintercept = 0, size = 1, alpha = 0.5)+
geom_line(aes(color = Condition, group = Condition))+
geom_point(size = 3.5, shape = 21, alpha = 0.5, aes(fill = Condition))+
theme_bw(base_size = 14)+
facet_grid(`Simulated Error`~., switch = "both")+
scale_y_continuous(limits = c(-3,3))+
xlab("CYP2D6 Genotype")+
ylab(expression(atop("Standardized Residuals", "("*CL[intrinsic]*")")))
CL_plot2.1CL_plot2.2 <- ggplot(data = sensitivity.CL,
aes(y = fct_reorder(Condition, standard_resid, .fun = median),
x = standard_resid))+
stat_boxplot(geom = "errorbar",
width = 0.5, size = 0.7) +
geom_boxplot(size = 0.7)+
geom_vline(xintercept = 0,
color = "blue", size = 1)+
geom_point(size = 2.5, shape = 21, aes(fill = Diplotype))+
# stat_summary(geom = "point",
# fun = mean,
# color = "red", size = 3)+
scale_x_continuous(limits = c(-3,3))+
theme_bw(base_size = 14)+
xlab(expression(atop("Standardized Residuals", "("*CL[intrinsic]*")")))+
ylab("Experimental Condition")+
geom_vline(xintercept = c(-2,2),
linetype = "dashed",
color = "darkgrey",
size = 1)+
ggeasy::easy_add_legend_title("CYP2D6\nGenotype")
CL_plot2.2plot(rich_covar.nlme, main = "Rich (0% CV)")plot(rich_CV5_covar.nlme, main = "Rich (5% CV)")plot(rich_CV10_covar.nlme, main = "Rich (10% CV)")plot(rich_CV20_covar.nlme, main = "Rich (20% CV)")plot(sparse3pt_covar.nlme, main = "Sparse (3pt, 0% CV)")plot(sparse3pt_CV5_covar.nlme, main = "Sparse (3pt, 5% CV)")plot(sparse3pt_CV10_covar.nlme, main = "Sparse (3pt, 10% CV)")plot(sparse3pt_CV20_covar.nlme, main = "Sparse (3pt, 20% CV)")plot(sparse4pt_covar.nlme, main = "Sparse (4pt, 0% CV)")plot(sparse4pt_CV5_covar.nlme, main = "Sparse (4pt, 5% CV)")plot(sparse4pt_CV10_covar.nlme, main = "Sparse (4pt, 10% CV)")plot(sparse4pt_CV20_covar.nlme, main = "Sparse (4pt, 20% CV)")plot(PM_covar.nlme, main = "PM Rich (0% CV)")plot(PM_CV5_covar.nlme, main = "PM Rich (5% CV)")plot(PM_CV10_covar.nlme, main = "PM Rich (10% CV)")plot(PM_CV20_covar.nlme, main = "PM Rich (20% CV)")plot(PM_3pt_covar.nlme, main = "PM Sparse (3pt, 0% CV)")plot(PM_3pt_CV5_covar.nlme, main = "PM Sparse (3pt, 5% CV)")plot(PM_3pt_CV10_covar.nlme, main = "PM Sparse (3pt, 10% CV)")plot(PM_3pt_CV20_covar.nlme, main = "PM Spase (3pt, 20% CV)")plot(PM_4pt_covar.nlme, main = "PM Sparse (4pt, 0% CV)")plot(PM_4pt_CV5_covar.nlme, main = "PM Sparse (4pt, 5% CV)")plot(PM_4pt_CV10_covar.nlme, main = "PM Sparse (4pt, 10% CV)")plot(PM_4pt_CV20_covar.nlme, main = "PM Sparse (4pt, 20% CV)")plot(rich_covar.nlme,log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (0% CV)",
xlab = "log(Fitted values)")plot(rich_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (5% CV)",
xlab = "log(Fitted values)")plot(rich_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (10% CV)",
xlab = "log(Fitted values)")plot(rich_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (20% CV)",
xlab = "log(Fitted values)")plot(sparse3pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 0% CV)",
xlab = "log(Fitted values)")plot(sparse3pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 5% CV)",
xlab = "log(Fitted values)")plot(sparse3pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 10% CV)",
xlab = "log(Fitted values)")plot(sparse3pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 20% CV)",
xlab = "log(Fitted values)")plot(sparse4pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 0% CV )",
xlab = "log(Fitted values)")plot(sparse4pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 5% CV)",
xlab = "log(Fitted values)")plot(sparse4pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 10% CV)",
xlab = "log(Fitted values)")plot(sparse4pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 20% CV)",
xlab = "log(Fitted values)")plot(PM_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (0% CV)")plot(PM_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (5% CV)")plot(PM_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (10% CV)")plot(PM_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (20% CV)")plot(PM_3pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt, 0% CV)")plot(PM_3pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt, 5% CV)")plot(PM_3pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt 10% CV)")plot(PM_3pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt, 20% CV)")plot(PM_4pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 0% CV)")plot(PM_4pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 5% CV)")plot(PM_4pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 10% CV)")plot(PM_4pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 20% CV)")An updated version of a previous function
population.sim() was generated called
population.sim2(). The updated function provides greater
flexibility substrate concentration range below 0.50 uM.
# ========== POPULATION SIMULATION FUNCTION =================#
## Updated to include low end point of 0.10 uM
population.sim2 <- function(
# Pre-define number of subjects per diplotype group
n = 10,
`CV%_D` = 25,
`CV%_V` = 0,
# Pre-define Michaelis Menten Kinetic Settings
Start = 0.5,
End = 2000,
By = 0.5,
# Turn off or on use of true population frequency
## If set to true "n" will equal total population number
pop_freq = FALSE,
# Pre-define reference data, and set seed
data = diplotype.data,
seed = 23457){
# ======== Create Data Containers ========
datasets <- list()
pop.data <- tibble()
# ======== Create Datasets by Diplotype ========
for (i in seq_along(data[[1]])){
# Setting Population Specific Frequency
n_pop <- if_else(pop_freq == FALSE, n,
ifelse(pop_freq == TRUE,
round(n*data[[i, c("Freq")]]), 0))
# Specify Diplotype for current loop
Diplotype = rep(data[[i,1]], n_pop)
# Random assignment of Vmax values N(0,sigma)
set.seed(seed)
Vmax_eta <- rnorm(n_pop, sd = `CV%_D`)
Vmax = data[[i,2]] + (Vmax_eta/100)*data[[i,2]]
# Random assignment of Km values N(0,sigma)
set.seed(seed)
Km_eta <- rnorm(n_pop, sd = `CV%_D`)
Km = data[[i,3]] + (Km_eta/100)*data[[i,3]]
# Store each diplotype group in a list container
temp <- tibble(Diplotype, Vmax, Km)
datasets[[i]] <- temp
}
# ======== Combine Diplotype Datasets ========
for (i in seq_along(data[[1]])){
pop.data <- rbind(pop.data, datasets[[i]])
}
# ======== Simulate Michaelis Menten ========
set.seed(seed)
pop.data <- pop.data %>%
mutate(ID = factor(seq_along(Diplotype),
levels = seq_along(Diplotype))) %>%
relocate(ID) %>%
expand_grid(S = c(0.1, seq(Start, End, By))) %>%
mutate(V = round((Vmax*S)/(Km+S), digits = 2))%>%
group_by(ID) %>%
# Additon of Residual Error
mutate(resid_pct = round(rnorm(n = length(S),sd = `CV%_V`),
digits = 3),
V_adj = V+(V*resid_pct/100)) %>%
ungroup()
return(pop.data)
}Michaelis-Menten kinetics were simulated for each subject in the virtual population across the 4 experimental designs (0%, 5%, 10%, and 20% CV).
# ========== CREATE BOOTSTRAP VIRTUAL POPULATIONS ============#
# Virtual Populations
CV0_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 0, End = 2000)
CV5_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 5, End = 2000)
CV10_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 10, End = 2000)
CV20_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 20, End = 2000)
# Sample Population Containers
CV0_Sample_populations = list()
CV5_Sample_populations = list()
CV10_Sample_populations = list()
CV20_Sample_populations = list()Using a for() loop, a total of 10 bootstrapped
populations were generated for each experiental design (0-20% CV) by
randomly sampling (with replacement) 1% of the corresponding virtual
population (10 subject per genotype group, n = 90).
## Bootstrap for-loop
B = 10 ## number of bootstraps
for (i in 1:B) {
# Bootstrap for 0% CV Population ---------
CV0_Sample_populations[[i]] <- CV0_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
# Bootstrap for 5% CV Population ---------
CV5_Sample_populations[[i]] <- CV5_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
# Bootstrap for 10% CV Population ---------
CV10_Sample_populations[[i]] <- CV10_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
# Bootstrap for 20% CV Population ---------
CV20_Sample_populations[[i]] <- CV20_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
}
# Combining Population Data -----------
Bootstrap.populations <- list(CV0_Sample_populations,
CV5_Sample_populations,
CV10_Sample_populations,
CV20_Sample_populations)
names(Bootstrap.populations)<- c("CV0_Pops",
"CV5_Pops",
"CV10_Pops",
"CV20_Pops")
# saveRDS(Bootstrap.populations, "R Output/Bootsrap Populations Data.rds")Strategic sampling designs were implemented using a custom function
update_pop(), which functions similarly to the
update_data() function mentioned previously, however, is
designed to work more efficiently the Bootstrap.populations
data which is a list of dataframes). The inputs for the
update_pop() function include the following:
pop_data - A dataframe or indexed dataframe from a list of dataframes.
type - Specifies the strategic sampling design to be implemented. Options include: “rich”, “sparse3pt”, “sparse4pt”, “PM”, “PM_3pt”, “PM_4pt”.
#======== CREATE BOOTSTRAP EXPERIMENTAL DATASETS =========#
# Function to transform virtual population into sample population
update_pop <-function(pop_data, type){
# Sampling Information --------------
# Full Range Set
full_set <- c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)
PM_range <- c(1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, 2000)
# Strategic Sampling Sets ----
sparse3pt_set1 <- c(0.1, 2.5, 50)
sparse3pt_set2 <- c(1, 10, 100)
sparse4pt_set1 <- c(0.1, 2.5, 25, 50)
sparse4pt_set2 <- c(1, 10, 50, 100)
PM_3pt_set1 <- c(10,100,1000)
PM_3pt_set2 <- c(25,250,2000)
PM_4pt_set1 <- c(1,10,100,1000)
PM_4pt_set2 <- c(5,25,250,2000)
PM_range_set1 <- c(1,10,25,50,100,400,1000,2000)
PM_range_set2 <- c(5,15,30,60,90,200,800,1600)
ID_key <- pop_data %>%
select(Diplotype, ID) %>%
unique() %>%
group_by(Diplotype) %>%
mutate(n = 1:n()) %>%
ungroup()
PM_Diplotype <- ID_key %>%
select(Diplotype) %>%
filter(Diplotype %in% c("4/5","4/41")) %>%
unique()
EM_Diplotype <- ID_key %>%
select(Diplotype) %>%
filter(!Diplotype %in% c("4/5","4/41")) %>%
unique()
# Rich Sampling Simulation ----------------
Data <- pop_data %>%
# Filter S Based on Experimental Model Type
filter(S %in% switch(type,
"rich" = full_set,
"sparse3pt" = full_set,
"sparse4pt" = full_set,
"PM" = PM_range,
"PM_3pt" = PM_range,
"PM_4pt" = PM_range)) %>%
#Assign Sampling Key
left_join(ID_key, by = c("Diplotype", "ID")) %>%
# Strategic Sampling
mutate(Set = if_else(n %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1",
S %in% switch(type,
"rich" = full_set,
"sparse3pt" = sparse3pt_set1,
"sparse4pt" = sparse4pt_set1,
"PM" = PM_range_set1,
"PM_3pt" = PM_3pt_set1,
"PM_4pt" = PM_4pt_set1),
S %in% switch(type,
"rich" = full_set,
"sparse3pt" = sparse3pt_set2,
"sparse4pt" = sparse4pt_set2,
"PM" = PM_range_set2,
"PM_3pt" = PM_3pt_set2,
"PM_4pt" = PM_4pt_set2))) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V, Set) %>%
# Refine Inidvidual in dataset according to model
filter(Diplotype %in% switch(type,
"rich" = EM_Diplotype[[1]],
"sparse3pt" = EM_Diplotype[[1]],
"sparse4pt" = EM_Diplotype[[1]],
"PM" = PM_Diplotype[[1]],
"PM_3pt" = PM_Diplotype[[1]],
"PM_4pt" = PM_Diplotype[[1]]))
# Data Output ----------
Data
}Using a for() loop, each virtual population was subject
to all 3 sampling strategies (rich, sparse 3pt, and sparse 4pt); with
sampling range varying according to CYP2D6 genotype as described
previously (see Strategic Sampling Approach). Data was combined
as a list of a list of dataframes with the first level (n = 4) being %CV
specific data, the second level (n = 3) being sampling strategy with
each variable at this level containing data for 10 unique population (n
= 120 total datasets).
# Create Experimental Data sets Populations ---------
CV0_Boot.data <- list()
CV5_Boot.data <- list()
CV10_Boot.data <- list()
CV20_Boot.data <- list()
# For Loop to apply strategic sampling to all dataframes ----------
for (i in 1:B) {
CV0_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "rich")
CV0_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "sparse3pt")
CV0_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "sparse4pt")
CV0_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "PM")
CV0_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "PM_3pt")
CV0_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "PM_4pt")
cat('CV0 Iteration', i, "of", B, "complete...\n")
CV5_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "rich")
CV5_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "sparse3pt")
CV5_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "sparse4pt")
CV5_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "PM")
CV5_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "PM_3pt")
CV5_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "PM_4pt")
cat('CV5 Iteration', i, "of", B, "complete...\n")
CV10_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "rich")
CV10_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "sparse3pt")
CV10_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "sparse4pt")
CV10_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "PM")
CV10_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "PM_3pt")
CV10_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "PM_4pt")
cat('CV10 Iteration', i, "of", B, "complete...\n")
CV20_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "rich")
CV20_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "sparse3pt")
CV20_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "sparse4pt")
CV20_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "PM")
CV20_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "PM_3pt")
CV20_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "PM_4pt")
cat('CV20 Iteration', i, "of", B, "complete...\n")
}
# Combine all boostrap data into a single data structure ---------
Complete_Bootstrap.data <- list(CV0_Boot.data,
CV5_Boot.data,
CV10_Boot.data,
CV20_Boot.data)
names(Complete_Bootstrap.data) <-c("CV0","CV5","CV10","CV20")
# Save the completed bootstrap data sets ------
# saveRDS(Complete_Bootstrap.data, "R Output/Complete Bootstrap Data.rds")fit_nlme() is a custom function used to simultaneously
fit non-linear least squares, and non-linear mixed effect models (with
and without covariates). The inputs for the function include the
following:
data - Data to be used for model fitting
type - Type of CYP2D6 metabolizer. Options include: “EM” for extensive metabolizer and ultra-rapid metabolizer data, and “PM” for intermediate and poor metabolizer data (default = “EM”).
which - Specifies the desired model output. Options include: “nls” = non-linear least squares, “base.model” for NLME with out covariates, “covar.model” for NLME with covariates (CYP2D6 genotype); default = “covar.model”.
# ==== Fitting Mixed Effect Model ================#
fit_nlme <- function(data, type = "EM", which = "covar.model"){
n <- if_else(type == "PM", 1, 6)
# Initial Screen ----
data0.grp <- groupedData(V~S|ID, data)
fit0.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = data0.grp)
# Extract Erroneous Subjects
error.vmax <- coef(fit0.nls) %>%
filter(is.na(Vmax)) %>%
row.names() %>%
as.numeric()
error.km <- coef(fit0.nls) %>%
filter(is.na(Km)) %>%
row.names() %>%
as.numeric()
error.list <- c(error.vmax, error.km)
# Clean Data set
data.grp <- data0.grp %>% filter(!ID %in% error.list)
# Base Model Fits
fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = data.grp)
base.nlme <- nlme(fit.nls, random = pdDiag(Vmax + Km ~ 1))
covar.nlme <- update(base.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(fixef(base.nlme)[[1]], rep(0,n),
fixef(base.nlme)[[2]], rep(0,n)),
weights = varPower())
# Output ----------
switch (which,
"nls" = fit.nls,
"base.model" = base.nlme,
"covar.model" = covar.nlme)
}extract_nlme_est() is a custom function used to fit,
extract, and tidy NLME parameter estimates from experimental data. This
function is optimized to work within the context of a for()
loop, and contains the fit_nlme() function as a dependency.
The inputs for the function include the following:
exp.data - Experimental data that NLME estimates are desired from.
condition - Label used to tag experimental condition from which the experimental data belongs.
PM - TRUE or FALSE argument specifying where data belongs to IM/PM CYP2D6 genotype (TRUE), or UM/EM CYP2D6 genotype (FALSE); default = “FALSE”.
pop.num - Label specifying which bootstrapped population is being analyzed (default = 1).
extract_nlme_est <- function(exp.data, condition, PM = FALSE, pop.num = 1){
# Input ----------
model.type <- if_else(PM == TRUE, "PM", "EM")
nlme.model <- fit_nlme(exp.data, type = model.type)
int.data <- intervals(nlme.model, which = "fixed")
label <- if_else(PM == TRUE, "4/41", "1/1")
# Tidy Data -------
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(Parameter = rownames(int.data$fixed)) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = paste0("Diplotype",label))) %>%
group_by(Parameter) %>%
mutate(est. = if_else(Diplotype != paste0("Diplotype",label),
est. + est.[Diplotype == paste0("Diplotype",label)],
est.),
across(where(is.numeric), round, 2)) %>%
ungroup() %>%
separate(Diplotype, remove = T, sep = "Diplotype", into = c(".", "Diplotype"))%>%
select(Parameter, Diplotype, est.)
colnames(temp) <- c("Parameter", "Diplotype", paste(condition))
# Output ----
temp %>%
pivot_longer( 3, names_to = "Condition", values_to = "Estimate") %>%
mutate(Population = pop.num)
}
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$sparse3pt[[3]],
condition = "Sparse (3pt, CV 5%)", PM = FALSE, pop.num = 3)Using a for() loop, bootstrapped parameter estimates
were extracted from experimental data
(Complete_Bootstrap.data) using the
extract_nlme.est() function. Parameter estimate for all
experimental designs and conditions for each population were stored in
object Complete_Bootstrap.est as a list. A collapsed
version was saved as Complete Bootstrap table.RDS for data
analysis and user reproducibility purposes.
Complete_Bootstrap.est <- list()
for (i in 1:B) {
# CV 0% Condition -------------
CV0_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$rich[[i]],
condition = "Rich (CV 0%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$sparse3pt[[i]],
condition = "Sparse (3pt, CV 0%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$sparse4pt[[i]],
condition = "Sparse (4pt, CV 0%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$PM[[i]],
condition = "Rich (CV 0%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$PM_3pt[[i]],
condition = "Sparse (3pt, CV 0%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$PM_4pt[[i]],
condition = "Sparse (4pt, CV 0%)", PM = TRUE, pop.num = i))
cat('CV0% - Iteration', i, "of", B, "complete...\n")
# # # CV 5% Condition -------------
CV5_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$rich[[i]],
condition = "Rich (CV 5%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$sparse3pt[[i]],
condition = "Sparse (3pt, CV 5%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$sparse4pt[[i]],
condition = "Sparse (4pt, CV 5%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$PM[[i]],
condition = "Rich (CV 5%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$PM_3pt[[i]],
condition = "Sparse (3pt, CV 5%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$PM_4pt[[i]],
condition = "Sparse (4pt, CV 5%)", PM = TRUE, pop.num = i))
cat('CV5% - Iteration', i, "of", B, "complete...\n")
# # CV 10% Condition -------------
CV10_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$rich[[i]],
condition = "Rich (CV 10%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$sparse3pt[[i]],
condition = "Sparse (3pt, CV 10%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$sparse4pt[[i]],
condition = "Sparse (4pt, CV 10%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$PM[[i]],
condition = "Rich (CV 10%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$PM_3pt[[i]],
condition = "Sparse (3pt, CV 10%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$PM_4pt[[i]],
condition = "Sparse (4pt, CV 10%)", PM = TRUE, pop.num = i))
cat('CV10% - Iteration', i, "of", B, "complete...\n")
# CV 20% Condition -------------
# Problematic fails to converge, requires further attention ##
CV20_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$rich[[i]],
condition = "Rich (CV 20%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$sparse3pt[[i]],
condition = "Sparse (3pt, CV 20%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$sparse4pt[[i]],
condition = "Sparse (4pt, CV 20%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$PM[[i]],
condition = "Rich (CV 20%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$PM_3pt[[i]],
condition = "Sparse (3pt, CV 20%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$PM_4pt[[i]],
condition = "Sparse (4pt, CV 20%)", PM = TRUE, pop.num = i))
cat('CV20% - Iteration', i, "of", B, "complete...\n")
Complete_Bootstrap.est[[i]] <- rbind(CV0_est.table,
CV5_est.table,
CV10_est.table,
CV20_est.table)
# Progress Report -------
cat('==== ROUND', i, "of", B, "COMPLETE! ===\n")
}
Complete_Bootstrap.table <- bind_rows(Complete_Bootstrap.est)
# saveRDS(R Output/Complete Bootstrap Table.rds")Non-parametric bootstrap analysis based on a total of 120 generated
data sets (10 re-sampled populations per experimental condition) was
performed to evaluate the precision of the final model parameters (see
separate Bootstrap Analysis section above). Mean parameter
estimates and 95% confidence intervals for each condition were
summarized by genotype. Confidence intervals were calculated using the a
custom confidence interval function CI() detailed below
based on (Equation 6); where confidence interval (CI) is equal to the
mean value of the sample population plus or minus the t-score at
significance level(0.05) for N-1 degrees of freedom multiplied by the
standard error of the sample population mean. The inputs for the
CI() function are as followed:
x - A vector of values for which a confidence interval is to be calculated.
alpha - Alpha level for confidence interval (default = 0.05)
which - Specifies the desired output. Options: confidence interval rang = “ci”, lower interval estimate = “lower”, upper interval estimate = “upper”, mean estimate = “est”, all estimates (mean [lower - upper]) = “all”.
# Import Bootstrap Data ----------
Bootstrap.table <- read_rds("R Output/Complete Bootstrap Table.rds")
# function to calculate confidence interval of bootstrap estimates -----
CI <- function(x, alpha = 0.05, which = "ci"){
mean <- mean(x)
n <- length(x)
std_dev <- sd(x)
std_error <- std_dev/sqrt(n)
degrees_of_freedom = n - 1
t_score <- qt(p=alpha/2, df=degrees_of_freedom, lower.tail=F)
margin_error <- t_score * std_error
lower <- mean - margin_error
upper <- mean + margin_error
ci <- paste0("(", round(lower, 2),
" - ",round(upper, 2),")")
all <- paste0(round(mean, 2),
" (", round(lower, 2),
" - ",round(upper, 2),")")
# Output ------
switch (which,
"ci" = ci,
"lower" = lower,
"upper" = upper,
"est" = mean,
"all" = all)
}Summary.boot <- Bootstrap.table %>%
mutate(Index = 1:length(Population)) %>%
filter(Index != 315) %>% # un-physiologic scenario
group_by(Parameter, Diplotype, Condition) %>%
summarize(Est = round(mean(Estimate), 2),
sd = round(sd(Estimate), 2),
Lower = CI(Estimate, which = "lower"),
Upper = CI(Estimate, which = "upper"),
`CI (95%)` = CI(Estimate, which = "ci"),
all = CI(Estimate, which = "all"),
n = length(Estimate)) %>%
ungroup() %>%
left_join(actual.est, by = c("Diplotype", "Parameter")) %>%
mutate(Residual = Actual - Est)# Vmax Plot --------------
ggplot(data = Summary.boot %>%
filter(Parameter == "Vmax"),
aes(x = Est, y = Condition))+
geom_vline(aes(xintercept = Actual),
color = "blue", size = 0.7, alpha = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.25)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*0.8)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.3)), alpha = 0)+
geom_vline(aes(xintercept = c(Actual*0.7)), alpha = 0)+
geom_errorbar(aes(xmin = Lower, xmax = Upper), width = 0.6)+
geom_pointrange(aes(xmin = Lower, xmax = Upper))+
facet_wrap(~Diplotype, scales = "free_x", ncol = 3, nrow = 3)+
theme_bw()+
ggtitle("Bootstrapped Confidence Intervals (95%) of Vmax")+
xlab(bquote(paste('V'['max']*" Estimate")))+
ylab("Experimental Condition")Bootstrapped parameter estimates with 95% confidence intervals (CI) for final model across all experimental conditions (stratified by CYP2D6 genotype). Black markers with intervals represent the mean estimate and 95% confidence interval of 10 bootstrapped sample populations, blue line represents true estimate of Vmax for specified genotype group, dotted red lines represent 80%-125% lower-and-upper bounds of parameter estimate.
# Km Plot ---------
ggplot(data = Summary.boot %>%
filter(Parameter == "Km"),
aes(x = Est, y = Condition))+
geom_vline(aes(xintercept = Actual),
color = "blue", size = 0.7, alpha = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.25)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*0.8)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.3)), alpha = 0)+
geom_vline(aes(xintercept = c(Actual*0.7)), alpha = 0)+
geom_errorbar(aes(xmin = Lower, xmax = Upper), width = 0.6)+
geom_pointrange(aes(xmin = Lower, xmax = Upper))+
facet_wrap(~Diplotype, scales = "free_x", ncol = 3, nrow = 3)+
theme_bw()+
ggtitle("Bootstrapped Confidence Intervals (95%) of Km")+
xlab(bquote(paste('K'['m']*" Estimate")))+
ylab("Experimental Condition")Bootstrapped parameter estimates with 95% confidence intervals (CI) for final model across all experimental conditions (stratified by CYP2D6 genotype). Black markers with intervals represent the mean estimate and 95% confidence interval of 10 bootstrapped sample populations, blue line represents true estimate of Km for specified genotype group, dotted red lines represent 80%-125% lower-and-upper bounds of parameter estimate
# Session Information ------
sessionInfo()## R version 4.2.0 (2022-04-22 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.3.4 cowplot_1.1.1 scales_1.2.0 ggpubr_0.4.0
## [5] ggeasy_0.1.3 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.9
## [9] purrr_0.3.4 readr_2.1.2 tidyr_1.2.0 tibble_3.1.8
## [13] ggplot2_3.3.6 tidyverse_1.3.1 nlme_3.1-157
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.2 lubridate_1.8.0 insight_0.17.1 webshot_0.5.4
## [5] httr_1.4.3 tools_4.2.0 backports_1.4.1 bslib_0.4.0
## [9] sjlabelled_1.2.0 utf8_1.2.2 R6_2.5.1 DBI_1.1.2
## [13] colorspace_2.0-3 withr_2.5.0 tidyselect_1.1.2 gridExtra_2.3
## [17] emmeans_1.7.4-1 compiler_4.2.0 performance_0.9.0 cli_3.3.0
## [21] rvest_1.0.2 xml2_1.3.3 labeling_0.4.2 bayestestR_0.12.1
## [25] sass_0.4.2 mvtnorm_1.1-3 systemfonts_1.0.4 digest_0.6.29
## [29] minqa_1.2.4 rmarkdown_2.14 svglite_2.1.0 pkgconfig_2.0.3
## [33] htmltools_0.5.2 lme4_1.1-29 dbplyr_2.2.0 fastmap_1.1.0
## [37] highr_0.9 rlang_1.0.4 readxl_1.4.0 rstudioapi_0.13
## [41] jquerylib_0.1.4 farver_2.1.1 generics_0.1.2 jsonlite_1.8.0
## [45] car_3.1-0 sjPlot_2.8.10 magrittr_2.0.3 Matrix_1.4-1
## [49] parameters_0.18.1 Rcpp_1.0.8.3 munsell_0.5.0 fansi_1.0.3
## [53] abind_1.4-5 lifecycle_1.0.1 stringi_1.7.6 yaml_2.3.5
## [57] carData_3.0-5 MASS_7.3-56 grid_4.2.0 sjmisc_2.8.9
## [61] crayon_1.5.1 lattice_0.20-45 splines_4.2.0 ggeffects_1.1.2
## [65] haven_2.5.0 sjstats_0.18.1 hms_1.1.1 knitr_1.39
## [69] pillar_1.8.0 boot_1.3-28 estimability_1.3 ggsignif_0.6.3
## [73] effectsize_0.7.0 codetools_0.2-18 reprex_2.0.1 glue_1.6.2
## [77] evaluate_0.15 modelr_0.1.8 nloptr_2.0.3 vctrs_0.4.1
## [81] tzdb_0.3.0 cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [85] datawizard_0.4.1 cachem_1.0.6 xfun_0.31 xtable_1.8-4
## [89] broom_0.8.0 coda_0.19-4 rstatix_0.7.0 viridisLite_0.4.0
## [93] ellipsis_0.3.2